Appendix A 

# 

# 

# File: trace.py 
# 

# Implements a ray tracing technique and related interpolation and 

# location procedures on 3D tetrahedral meshes. This class is 

# closely related to the Geometry class. This example is based on 

# a tet mesh, but the same functions could be provided for other 

# mesh types. 
# 

# Three major example functions are provided: 

# 1) cellNunib - Given a set of (x,y,z) coordinates, finds the cell 

# which contains those coordinates. This could be used for 

# example to assign a CT voxel with known coordinates 

# to a cell on the tet mesh. 

# 2) fieldVal - Given a set of (x,y,z) coordinates and field values 
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# at the mesh cell vertexes, determines the LDFEM field value 

# at those coordinates. This could be used in a post 

# processing utility to provide exact values at specified 

# points in the problem. 

# 3) trackData - Given start and end (x,y,z) coordinates of a 

# track, makes a list of the mesh cells and the path lengths in 

# each cell that the the track passes through. This could be 

# useful in a euialytic ray trace algorithm to calculate 

# total optical path length along a ray when cells have 

# different materials assigned to them. The ray starting 

# point must lie inside the mesh, but the ending point can 

# be external. This is useful when modeling external point 

# sources for example. 
# 

# Searches through the mesh are expedited in this class by 

# employing a line search based algorithm that reduces the number 

# of cells to be searched. This algorithm effectively bins the 

# mesh cells. Other binning-like algorithms could be envisioned. 

# A production algorithm would include numerous tests and error 

# traps that are not included herein for the sake of simplicity. 
# 

# 



# Author: John McGhee, Radion Technologies 
# 

# Date : 09 March 2004 

# 

# 

# $Id: trace. py,v 1.3 2004/03/10 17:17:26 mcghee Exp $ 
# 

# 

from sys import exit 
from math import sqrt 

class TetTrace: 

"Data and methods related to tracing, interpolation, euid point location on 3D 

tet meshes" 

def init (self ,meshdata) : 

"Initializes the TetTrace class" 

self .meshdata = meshdata 
self. path = [] 
self . end_vol_coord = [] 
self . end_t race = 0 
self . start_coord = [] 
self . end_coord = [] 
self. omega = [] 
self . total_path = 0. 

# 

def f ieldVal {self , xyz_coord, field_data) : 

"Finds the cell containing an arbitrary mesh point, and " + \ 

"returns the LDFEM field value at that point" 
icell = 0 

end_coord = xyz_coord 

start_coord = self . centroid( icell) 

i = self .trackData (icell, start_coord, end_coord) 

return self . interp (f ield_data) 

# 

def interp(self, field_data) : 

"Interpolates to find LDFEM field value inside a tet" 
icell = self .patht-U [0] 
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xO = field_data[icelll [0] 
xl = field__data[icell] [1] 
x2 = field_data[icell] [2] 
x3 = field_data[icell] [3] 

return ( xO*self . end_vol_coord[0] + xl*self . end_vol_coord[l] + 

x2*self .end_vol_coord[2] + x3*self . end_vol_coord[3] ) 



def cellNuinb(self , xyz_coord) : 

"Finds the cell containing an arbitrary mesh point" 
icell = 0 

end_coord = xyz_coord 

start_coord = self . centroid( icell) 

i = self . trackData (icell, start_coord, end_coord) 

return self .path (-1] [0] 

# 

def trackData (self , start__cell, start_coord, end_coord) : 

"Calculates the cell nxambers and path lengths in a ray" 

self. path = [] 

self . end_cell = -999 

self . end_vol_coord = [] 

self . end_trace = 0 

self . star t_coord = start_coord[ : ] 

self . end__coord = end_coord[:] r 

> 

self. omega = [ end_coord[0] -start_coord[0] , 
end_coord[l] -start_coord[l] , 
end_coord[2] -start_coord[2] ] 
self . total_path = sqrt (self .omega [0] **2 + self .omega [1] **2 + 

self .omega [2] **2) 
if (self . total_path > 0.): 
XX = 1 . /self . total_j)ath 

self .omega [0]=xx*self. omega [0] ; self .omega [1] =xx*self. omega [1] 
self. omega [ 2 ] =xx* sel f , omega [ 2 ] 

i count = -1 

exit_coord = self. start_coord[ : 1 
next_cell = start_cell 

# walk through the ray one tet at a time, 
while (not self . end^trace) : 

icount = icount+1 

if {icount >= self .meshdata.ncells-1) : 

print "Error: cells in track exceeds total cells in mesh" 

exit ( ) 

entry_coord=exit_coord [ : ] 
exit_coord = [ ] 
current_cell = next_cell 

next_cell, exit_coord = self . cell_path {current_cell , entry_coord) 

# print self .path [-1] #list the cell numbers and path lengths for debug 

# Perform some error checking. 
X = 0. 

for i in self .path: 

if (i[0] < 0) or (i[03 >= self .meshdata.ncells) : 

print "Error: cell index out of range in trackData" 
exit ( ) 
X = x + i [1] 
if (abs(l. - x/self .total_path ) > l.e-6 ): 

print "Error: total path length in error in trackData" 
exit () 
return len( self .path) 
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def cell_path ( self , current_cell, entry__coord) : 

"Returns exiting coordinates and adjacent cell of a track through the 
current_cell " 

5 # ** Note that omega must have unit magnitude. ** 

# Check to see if the trace has exited the geometry 
if (current_cell == self .meshdata.bdry_f lag[0) ) : 

self . end__trace = 1 #true 
10 return (-999, (0., 0., 0.) ) 

# Check to see if the track terminates inside this cell. 

# Note that for computational efficiency if all the rays 

# terminate at one point, once the cell for that point is 
15 # found there is no need to perform repeated searches here. 

# A production algorithm would store this information 

# and avoid this test on every cell if the geometry permitted, 
vc = self .vol_coord( self . end_coord, current_cell) 

incell = 1 #true 
20 for i in vc : 

if (i < 0. ) or (i > 1. ) : 
incell = 0 #false 
break 
if (incell) : 

25 cell_path = sqrt( (entry_coord[01 -self . end_coord[ 0] ) **2 + \ 

(entry_coord[ll-self .end^coordCl] ) **2 + \ 
(entry_coord[2]-self .end_coord[2] ) **2 ) 
self .path. append ( (current_cell, cell_path) ) 
self .end_vol_coord = vc 
30 self .end_trace = 1 

return (-999, (0., 0., 0.) ) 

else : 

# Find the path length in this cell and the exiting face. 
35 cell_path = l.el5 

outface = -999 

for i in range (self .meshdata.nf aces) : 
xomp = [ ] 

iO = self .meshdata . f ace_node [i] [ 0] 
40 nO = self .meshdata. cell_node [cur rent_cell] [iO] 

xO = self .meshdata. node_coord[nO] [0] 
yO = self .meshdata .node_coord[nO] [1] 
zO = self .meshdata .node_coord[nO] [2] 
xomp . append ( xO - entry_coord [ 0 ] ) 
45 xomp . append (yO - entry_coord [ 1 ] ) 

xomp . append ( z 0 - entry_coord [ 2 ] ) 
avect = self .meshdata. avect (current_cell, i) 
omega_dot_a = avect [0] *self . omega [0] + \ 

avect [1] *self .omega [1] + avect [2] *self . omega [2] 
50 if (omega_dot_a > -l.e-30): 

omega_dot_a = max ( omega_dot_a , l.e-30) 
xomp_dot_a = avect [0] *xomp[0] + \ 

avect [1] *xomp[l] + avect [2 ] *xomp [2 J 
if (xomp_dot_a > -l.e-30): 
55 xomp_dot_a = max (xomp_dot_a, l.e-30) 

XX = xomp_dot_a/omega_dot_a 
if (XX < cell_path) : 
outface = i 
cell_path = XX 
60 if (outface < 0) : 

print "Error: unable to find an exiting face" 
exitO 

exit_coord = ( entry_coord [ 0 ] + cell_path* sel f . omega [ 0 ] , 
entry_coord [1] + eel l_path* sel f . omega [1] , 
65 entry^coord [2 ] + cell_path*self ,omega[2 3 ) 

self .path. append ( (current_cell, cell_path) ) 
next_cell = self .meshdata. adjcell (current^cell, outface) [0] 
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return (next_cell, exit__coord) 



de£ centroid(sel£, icell) : 

"Returns the centroid of cell icell" 



nO 




self. 


. meshdata , 


.cell. 


.node [icell] [0] 


nl 




sel f . meshdata , 


.cell. 


.node [icell] [1] 


n2 




self, 


.meshdata. 


.cell_node[ icell] [2] 


n3 




self, 


. meshdata . 


. cell. 


.node [icell] [3] 


xO 




sel f . meshdata . 


. node. 


.coord [nO] [0] 


yo 




self, 


. meshdata , 


. node. 


.coord [nO] [1] 


zO 




self .meshdata . 


. node. 


.coord [nO] [2] 


xl 




self, 


.meshdata , 


. node. 


.coord [nl] [0] 


yi 




self, 


. meshdata , 


. node. 


_coord[nl] [1] 


zl 




self. 


.meshdata , 


. node. 


_coord[nl] [2] 


x2 




self. 


.meshdata. 


. node. 


.coord [n2] [0] 


y2 




self. 


.meshdata , 


. node. 


.coord [n2] [1] 


z2 




sel f . meshdata . node. 


.coord [n2] [2] 


x3 




self 


. meshdata . 


. node. 


.coord [n3] [0] 


y3 




self, 


. meshdata . 


, node. 


.coord [n3] [1] 


z3 




self. 


•meshdata, 


. node. 


coord [n3] [2] 



return ( 0 . 25* (x0+xl+x2+x3 ) , 
0.25* (y0+yl+y2+y3) , 
0.25*(z0+zl+z2+z3) ) 



def vol_coord(self , xyz_coord/ icell) : 

"Returns the icell based volume coordinates of a point" 



nO 




self. 


.meshdata . 


cell_ 


node [icell] [0] 


nl 




self. 


. meshdata . 


cell_ 


.node [icell] [1] 


n2 




self. 


. meshdata . 


cell.. 


node [icell] [2] 


n3 




self. 


. meshdata . 


cell_ 


.node [icell] [3 ] 


xO 




self, 


. meshdata . 


node_ 


.coord [nO] [0] 


yo 




self. 


. meshdata . 


node_ 


.coord [nO] [1] 


zO 




self. 


.meshdata . 


node_ 


.coord [nO] [2] 


xl 




self. 


.meshdata. 


node_ 


.coord [nl] [0] 


yi 




self .meshdata . 


node_ 


coord [nl] [1] 


zl 




self . meshdata . 


node_ 


.coord[nl] [2] 


x2 




self. 


. meshdata . 


node_ 


.coord [n2] [0] 


y2 




self . meshdata . 


node_ 


.coord [n2] [1] 


z2 




self. 


. meshdata . 


node_ 


.coord [n2] [2] 


x3 




self. 


. meshdata . 


node_ 


.coord [n3] [0] 


y3 




self 


. meshdata . 


node_ 


.coord [n3] [1] 


z3 




self. 


.meshdata . 


node_ 


.coord [n3] [2] 


amat 


= [ 


[ n, [], 


[] ] 










[ [], [], 


[] ] 










[ [], [], 


[] ] 


] 



amat [0] [0] 




Xl 


- xO 


amat [0] [1] 




x2 


- xO 


amat [0] [2] 




x3 


- xO 


amat [1] [0] 




yi 


- yO 


amat[l] [1] 




y2 


- yO 


amat [1] [2] 




y3 


- yO 


amat [2] [0] 




zl 


- zO 


amat [2] [1] 




z2 


- zO 


amat [2] [2] 




z3 


- zO 


bvec = [ 0 . 


■ r 


0., 


0. ] 



bvec [ 0 ] = xyz_coord ( 0 ] - xO 
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bvec[l] = xyz_coord[l] - yO 
bvec[2] = xyz_coord[2] - zO 

self .ge3 {amat,bvec) 

XX = [] 

xx.appendd. - bvec[0] - bvec[l] -bvec[2]) 
XX . append ( bvec [ 0 ] ) 
XX . append ( bvec [ 1 ] ) 
XX . append ( bvec [ 2 ] ) 

return xx 



ge3(self, amat, bvec): 

"Solves a 3x3 linear system using Grout's method with partial pivoting" 

# This is a general procedure that would likely 

# be obtained from a optimized linear algebra library in a 

# production code. 

# j = 0 
imax = 0 

colmx = abs (amat [0] [0] ) 

dvmt = abs ( amat [ 1 ] [ 0 ] ) 
if (dum >= colmx) : 

imcLx = 1 

colmx = dum 

dum = abs (amat [2] [0] ) 
if (dum >= colmx) : 

imax = 2 

colmx = dum 

if (imax != 0) : 

d\am = amat [imax] [0] 

amat [imax] [0] = amat[0][0] 

amat [ 0 ] [ 0 ] = dum 

dum = amat [imax] [1] 

amat [imax] [1] = amat[0][l] 

amat [ 0 ] [ 1 ] = dum 

dum = amat [imax] [2] 

amat [imax] [2] = amat[0][2] 

amat [ 0 ] [ 2 ] = dum 

dum = bvec [imax] 

bvec [ imax ] = bvec [ 0 ] 

bvec [ 0 ] = dum 

if (amat[0] [0] == 0.) : 

print "Error: zero pivot in ge3" 
exit () 

amat[0][0] = 1 . /amat [ 0] [0] 
amat[l](0] = amat [1] [0] *ainat [0] [0] 
amat[2][0] = amat [2] [0] *amat [0] [0] 

# j = 1 

amat[l][l] = amat[l][l] - amat [1] [0] *amat [0] [1] 
amat[2][l] = amat[2][l] - amat [2] [0] *amat [0] [1] 

imax = 1 

colmx = abs (amat [1] [1] ) 
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duin = abs (amat [2 ] [ 1] ) 
if (dum >= colmx) : 

imax = 2 

colmx = dum 

if (imax ! = 1) : 

dum = amat [imax] [0] 

amat[im£ix] [0] = amat[l] [0] 

amat [1] [0] = dum 

dum = amat [imax] [1] 

amat [imax] [1] = amat[l][l] 

amat[l] [1] = dum 

dum = amat [imax] [2 ] 

amat [imax] [2] - amat[l] [2] 

amat [1] [2] = dum 

dum = bvec[imax] 

bvec[imeLx] = bvec[l] 

bvec[l] = dum 

if (amat [1] [1] == 0. ) : 

print "Error: zero pivot in ge3" 
exit ( ) 

amat[l][l] = 1 . /amat [1] [1] 
amat[2][l] = amat [2] [1] *amat [1] [1] 

# j = 2 

amat[l][2] = amat[l][2] - amat [1] [0] *amat [ 0] [2] 

amat[2][2] = amat[2][2] - amat [2] [0] *amat [0] [2] - amat [2] [1] *amat [1] [2] 

if (amat [2] [2] == 0. ) : 

print "Error: zero pivot in ge3" 

exit ( ) 

amat[2][2] = 1 . /amat [2] [2] 

# forward cind backward substitution 
bvec[l] = bvec [1] -amat [1] [0] *bvec [0] 

bvec[2] = ( bvec[2]-amat[2] [0]*bvec[0]-amat[2] [l]*bvec[l] 
>*amat[2] [2] 

bvec[l] = ( bvec[l] -amat [1] [2] *bvec [2] ) *amat [1] [1] 

bvec[0] = ( bvec[0]-amat[0] [l]*bvec[l]-amat[0] [2]*bvec[2] )*amat[0] [0] 

return 1 

# 

# end of trace. py 

# 

# Sn Angular Quadrature Data 

# Set Name - Triangle Cliel>ycliev Lobatto, 18 angles 

sph__degree 
1 

end_sph_degree 

# cung# mu eta zi wgt 
quadrature 

1 -l.OOOOOOOE+00 O.OOOOOOOE+00 0 . OOOOOOOE+00 8 . 3333333E-02 

2 -4.4721360E-01 -8 . 2634298E-01 -3 . 4228247E-01 5 . 2083333E-02 

3 -4.4721360E-01 -3 . 4228247E-01 -8 . 2634298E-01 5 . 2083333E-02 

4 4.4721360E-01 -8 . 2634298E-01 -3 . 4228247E- 01^ 5 . 2083333E-02 

5 4.4721360E-01 -3 . 4228247E-01 -8 . 2634298E-01 5 . 2083333E-02 

6 -4.4721360E-01 8 . 2634298E-01 -3 . 4228247E-01 5 . 2083333E-02 
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encL.quadr a tur e 



# 

# 

# File: aquad.py 
# 

# Angular quadrature related methods and data. Any angular 

# quadrature can be used. These are the "ordinates" in the 

# discrete ordinates (Sn) method of angular discretization. 
# 

# 

# Author: John McGhee, Radion Technologies 
# 

# Date : 19 Feb 2004 

# 

# 

# $Id: aquad.py, V 1.7 2004/03/05 17:19:58 mcghee Exp $ 
# 

# 

from sys insert exit 

class Quadrature: 

"Provides information on the angular quadrature." 

def init (self , quad_f ile) : 

"Reads in the eungular quadrature from disk. " 



from string import split, atof, atoi, strip 
from sph import sphfun 

# Setup for read of geometry data from disk, 
se 1 f . name=quad_f i 1 e 

f = open (quad_f ile) 
record = f.readlineO 

# Read in the desired spherical harmonics degree: 
while (record and strip (record) != " sph_degr ee " ): 

record = f.readlineO 
if (record == " " ) : 

print "Error: unable to find sph_degree keyword on quadrature file 
%s\n" % quad_file 

exit () 
record = f . readline ( ) 
self . sph_degree = atoi (record) 
if (self . sph_degree < 1): 

print "Error: spherical harmonics degree (Imax) must be >= 1" 

exit () 

# Read in the quadrature points from disk. 

while (record and strip (record) != "quadrature" ) : 
record = f.readline() 
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if (record == " " ) : 

print "Error: unable to find quadrature keyword on quadrature file 
%s\n" % quad_file 

exit ( ) 
self .quad = [] 
record = f.readlineO 

while (record and strip (record) != "end^quadrature" ) : 
j = split (record) 
if (len(j) != 5): 

print "Error: wrong number of entries on quadrature data line" 

exitO 

self .quad. append (map (atof, j [1: 5] ) ) 
record = f . readline ( ) 
if (record ==""): 

print "Error: unable to find end_quadrature keyword on quadrature 
file %s\n" % quad_file 
exit ( ) 

self.nang = len( self .quad) 

# indexes between from the moment ordinal index and the (l,in) index 
self.nmom = (self . sph_degree+l ) **2 
self . index_!Lm = [] 
self .index_kk = {} 
k = -1 

for i in range ( self . sph_degree+l ) : 
for j in range ( -i, i+1) : 
k = k+1 

X = ( i , j ) 

self . index_lm. append (x) 
self . index_kk [x] = k 



# This demo uses the "standard" method 

# of conversion from discrete to moment and visa versa. 

# A "Galerkin" or other more advanced mappping 

# method may be useful for certain problem types with highly 

# forward peaked scattering. These advanced methods 

# would be implemented here. 



# discrete to moment array, 
self .dis2moin_data = [] 

for imom in range (self .nmom) : 
X = [] 

1 = self .index_lm[ imom] [0] 
m = sel f . index_lm [ imom] [1] 
for kk in range (self ,nang) : 
omega = sel f . qdpt ( kk ) 

X, append ( omega [3 ] *sphfun(l, m, omega [0], omega [1], 

omega [ 2 ] ) ) 
self .dis2mom_data. append (x) 

# moment to discrete array, 
self .mom2dis_data = [] 

for kk in range ( self . nang) : 
omega = self .qdpt (kk) 
X = [] 

for imom in range (self .nmom) : 
1 = self . index_lm [ imom] [ 0 ] 
m = self .index_lm[ imom] [1] 

X. append ( sphfund, m, omega [0], omega [1], omega [2 ] ) ) 
self .mom2dis_data . append (x) 

# 



def qdpt (self ,kk) : 

"Returns quadrature data for a single angle" 
if ( kk < 0 ) or ( kk >= self.nang ) : 

print "Error: cuigle index out of range." 
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exit () 
return self . quad [kk] 

# 

def dis2moin(self , kk) : 

"Returns a column of the discrete to moment array" 
if (kk < 0) or (kk > self.nang): 

print "Error: kk index out of range in dis2mom" 

exit ( ) 
X = [] 

for imom in range ( self . nmom) : 

x . append ( sel f . dis2mom_dat a [ imom] [ kk] ) 
return x 



# 

def mom2dis (self / kk) : 

"Returns a row of tlie moment to discrete array" 
if (kk < 0) or (kk > self.nang): 

print "Error: kk index out of range in mom2dis'' 

exit ( ) 

return self .mom2dis_data [kk] 

# 

# end of aguad.py 

# 



$Id: background_info. asc, V 1.4 2004/03/05 21:51:31 mcghee Exp $ 



1 . Introduction 



Tliis is an overview of a demonstration code to be included with Radion Technologies 
utility patent for the application of deterministic transport methods to medical 
therapy and imaging problems, 19 Mar 2004. The purpose of this code is to provide an 
example of the "preferred embodiment" of a 3D linearized Boltzmann trsuisport (BTE) 
solver, and to demonstrate that the algoritlim has been "reduced to practice". The 
demonstration code seeks to provide enough detail so that a worker "skilled in the 
art" could reproduce the significant aspects of the algorithm without difficulty. 
However, the demonstration code is not a complete production ready product, the 
scope is limited to demonstrating the key features involved. 



2 . Development of Mathematical Model 



Particle transport is a physical process. This physical 

process can be mathematically modeled l>y a integro-partial differential equation 
(PDE) . 

In order to construct this model we make a series of 

more-or-less standard assumptions regarding the physical process of particle 
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transport. These assumptions are as follows: 

particles are treated as points, no quantum mecheuiical wave effects are 
included. No wave effects such as polarization, refraction, or interference are 
included. 

particles travel in straight lines between point collisions, particles are 
unaffected by external forces such as gravity or electrical and magnetic fields. 

particle-particle interactions are neglected. The particle density 
is assumed to be low compared with the density of the surrounding medium, 
collisions are assumed to be instantaneous, scattering events are 
not delayed. 

-- material properties are isotropic, there is no preferred direction 
of particle travel within a single material. 

— material properties are known a priori and are not a function of 
the particle distribution. 

— it is assiimed that the particle density is sufficiently high that a mean value 
of the distribution adequately describes the solution. 

Under these assun^tions we construct what we refer to as the linearized Boltzmann 
Transport Equation (BTE) . This equation is a hyperbolic 
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integro partial differential equation. The solution to this equation provides a 
distribution of particles in space, angle, and energy which 

is an (hopefully accurate) approximation to the actual particle distribution 

associated with the physical transport process we are trying to model. 

For charged particle problems a continuous-slowing-down (CSD) term in energy from the 
Boltzmann-Fokker- Planck (BFP) equation is included in the BTE. 

We assume that only the steady state solution is of interest. The 

BTE has a term included for time-dependent problems, but we elect not 

to include this term in our formulation as the radiation equilibration times 

are very short compared to the time scales of interest for the present 

application. We have included this term on other projects without difficulty. 

3 . Development of Numerical Model 

Given the mathematical model described above, a discrete approximation to that 
model is created which is amenable to numerical solution. 

The algorithm we employ consists of a set of individual techniques extant in the 
transport literature for some time. However, these techniques are assembled into 
a novel and particularly effective combination for our application. The 
techniques used are: 

multi-group discretization in energy, with a discontinuous differencing in 
energy of the BFP CSD term. 

discontinuous finite element (DFEM) spatial differencing on an unstructured 
finite element mesh, and, 

— discrete ordinates (Sn) differencing in angle. 

4. Defining (critical) features of demonstration code. 
Equation solved: 

-- a discretized approximation to the first order form of the 

linearized BTE with an additional BFP CSD term for charged particles, a first 
order hyperbolic integro partial differential equation 



Discretization : 

— multi -group energy differencing 

— unstructured finite element mesh 
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— discontinuous finite element spatial differencing (DFEM) 

— discrete ordinates angular differencing (Sn) 

— spherical harmonics expansion of scattering source 

-- produces a large, sparse, asymmetric linear system of equations 

in a six dimensional phase space. 
Solution technique: 

Source iteration (Richardson iteration) 

— Sweep ordering and cycle breaking to produce lower triangular system -- 
Diffusion Synthetic Acceleration (DSA) of the source iteration 

process . 

Important features to be demonstrated: 

— multiple particle problems are supported 
interpolation of solution on FE mesh. 

— ability to restart from previous solution 

— ability to control refinement in space, angle, and energy 

Other significant novel aspects: 

— 3D lobatto-chebychev angular quadratures 

— ray tracing option to minimize ray effects for point sources, including 

first, second, and subsequent scattered sources. 

— memory management techniques 

— ability to provide adjoint solutions 

— ability to accommodate a variety of material properties databases in a transparent and 

modular fashion . 

— alternative solution methodologies - Krylov solvers. 

parallel ization: angular domain decomposition, sweep scheduling, 
Krylov solvers 

5 . Key Advantages of the Method Described 

Many techniques for discretizing the BTE have been developed over the 

years. The method described herein is an attractive solution algorithm because it 
provides a capability to capture the true transport solution everywhere in a problem with 
a minimal amount of computational effort. This method 

provides "the accuracy heretofore attributed exclusively to Monte Carlo, without the 
prohibitive computational cost." The major features which provide this capability are as 
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follows : 

— the Sn, multi-group method with arbitrary order anisotropic scatter is a convergent 
method. That is to say that it captures all the physics in 

the BTE and converges to the true transport solution as the computational mesh is refined in 
angle, space, and energy. 

arbitrary finite element mesh accurately captures the true problem 
geometry with a minimum number of spatial elements 
-- discontinuous finite element differencing captures 

discontinuities in solution and material properties and is 3rd order accurate, 
acceleratable, damped, and has the diffusion limit. 

the method provides a complete solution everywhere in the problem, free from 
statistical noise. There is no need to guess areas of interest 
beforehand as with Monte Carlo techniques . 

— rapid solutions are possible, arbitrary mesh provides capability to add 
resolution where needed and yet remain coarse elsewhere. DFEM spatial differencing 
remains accurate even on coarse meshes. 

— iterative solution technique provides a means to utilize information from 
previously solved problems to minimize cost of subsequent perturbation analysis. 

ability to control resolution in space, angle, and energy on a 
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problem by problem basis saves computational time by avoiding excess resolution in 
problems or problem regions which do not rec[uire it. 

6 . References 

1) Wareing, T.A. , J.M. McGhee, J.E. Morel, S.D. Pautz, "Discontinuous Finite 
Element Sn Methods on Three Dimensional Unstructured Grids", Nuclear Science and 
Engineering, Vol 138, 256-268 (2001) . 

2) Pautz, S.D., "An Algorithm for Parallel Sn Sweeps on Unstructured Meshes", 
Nuclear Science and Engineering, Vol 140, 111-136, (2002) . 

3) Lewis, E.E. and W.F. Miller, Jr. "Computational Methods of Neutron Transport", 
John Wiley and Sons, New York, 1984. 

4) Zienkiewicz, O.C. , and R. L. Taylor", "The Finite Element Method", Fourth 
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# 

# 

# File: bte.py 



# 

# Sets up and solves the Boltzmann transport equation for 

# a specific group, angle, and cell. This class is specific 

# to tets for simplicity of exposition. In a general purpose code 

# ' this class would likely be a virtual base class, with derived 

# classes for several different element types. 
# 

# 



# Author: John McGhee, Radion Technologies 
# 

# Date : 19 Feb 2004 

# 

# 

# $Id: bte.py, v 1.17 2004/03/05 17:19:58 mcghee Exp $ 
# 

# 

from sys import exit 

class BteEquation: 

"Forms the LD linear tet BoltzmcUin Transport Equation" 
cell_type = "tetrahedra" 

def init (self, geom, icell, icp, omega, sigt, beta, kk) : 

"Constructs the left hand side of the LD linear tet Boltzmann Equation" 

self .ib= [] 
self.bvec = [] 
self. icell = icell 
self. icp = icp 
self.n = 4 
self.qcsd = [] 
self . sum_qcsd = 0. 
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10 



15 



20 



25 



30 



self.gainma = 0. 
self.psiE = 0. 
self.kk = kk 

# Form omega dot A for each face in the cell 
X = (] 

for iface in range (geom. nf aces) : 
XX = 0 . 

yy = geom. avect (icell, if ace) 
for j in range (geom. ndim) : 

XX = XX + omega [ j ] *yy ( j 1 
XX = xx/12. 
X. append ( XX ) 

a = [] 
b = [] 
c = [] 
for i in x: 

if ( i < -geom. dplimit ) : 

a . append ( 0 . ) 

b . append ( abs { i ) ) 

c . append ( 1 ) 

elif ( i > geom. dplimit) : 

a . append ( i ) 

b . append ( 0 . ) 

c . append ( 0 ) 
else: 

a . append < 0 . ) 

b. append(0 . ) 

c . append ( 0 ) 



35 



40 



45 



50 



# Construct the removal term 
r = sigt*geom.vol ( icell) /20 . 

# Construct the CSD terms for charged particle problems. 

# Charged particle problems introduce an additional unlcnown 

# into the solution trial space, which introduces an 

# additional column and row into the transport matrix. We 

# elect eliminate the extra unknown from the matrix by cyclic 

# reduction, thereby keeping the number of rows and columns in 

# the transport matrix the same for both charged and neutral 

# particle problem types, 
if (icp) : 

self.kl = beta*geom, vol (icell) /4 . 

self.k2 = l./( 4.*(a[0]+a[l]+a[21+a[3] ) + \ 

sigt*geom. vol (icell) /3 . + 4.*self.kl) 
XX = self .kl*self .kl*self .k2 
cpO = XX + 0.2*self.kl 
cpl = XX + 0.4*self.kl 



else: 



CpO = 
cpl = 
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# Form the transport 


matrix 


for this cell 






self 


. amat= [ 


[] 


, [] ] 










self 


. amat [0] 


, append ( 


x[0] 


+ 


2.*(a[l]+a[2]+a[3] ) 


+ 


self 


. amat [0] 


. append ( 


x[0] 


+ 


a[2]+a[3] + r 


+ cpO 


) 


self 


. amat [ 0 ] 


. append ( 


x[0] 




a[l]+a[3] + r 


+ cpO 


) 


self 


. amat [0] 


. append ( 


x[0] 


+ 


a[l]+a[2] + r 


+ CpO 


) 


self 


. amat [ 1 ] 


. append ( 


x[13 


+ 


a(2]+a[3] + r 


+ cpO 


) 


self 


.amat [1] 


. append ( 


x[l] 


+ 


2.*(a[01+a[21+a[3] ) 


+ 


self 


.amat [13 


. append ( 


x[13 


+ 


a[0]+a[31 + r 


+ cpO 


) 


self 


.amat [1] 


. append ( 


x[l] 


+ 


a[0]+a[2} + r 


+ cpO 


) 


self 


.amat [2] 


. append { 


x[2] 


+ 


a[l]+a[3] + r 


+ cpO 


) 


self 


. amat [ 2 ] 


. append ( 


x[2] 


+ 


a[03+a[3] + r 


+ cpO 


) 


self 


. amat [ 2 ] 


. append ( 


x[2] 




2.*(a[0]+a[l]+a[3]) 


+ 



+ 2 , *r + cpl) 



2.*r + cpl) 



•'r + cpl) 
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self .amat [2] .append ( x[2] + a[03+a[l] + r + cpO ) 

self .amat [3] .append ( x[3] + a[l]+a[2] + r + cpO ) 

self , amat [3] .append ( x[3] +a[0]+a[2] +r+ cpO ) 

self . amat [3] . append( x[3] + a[03+a[l] + r + cpO ) 

self .amat [3] .append ( x[3] + 2 . * (a[0]+a[l] +a[2] ) + 2.*r + cpl) 

# Construct the coefficients used to form the in-flow terms, 
self .yy= [] 

self . yy . append ( 

( 2.*b[l], 2.*b[23, 2.*b[3], 0., b[2], b[3] , 
0., btl], b[3], 0., b[l], b[21 ) ) 
self. yy . append ( 

( 0., b[2], b[3], 2.*b[0], 2.*b[2], 2.*b[3], 
b[0]. 0., b[3], b[0], 0., b[2] ) ) 
self . yy . append ( 

( btl], 0., b[3], b[0], 0., b[31, 2.*b[0], 
2.*bCll, 2.*b[3], b[0], b[ll, 0.) ) 
self. yy . append ( 

( b[l], b[2], 0., b[0], b[2], 0., b[0], b[ll , 
0., 2.*b[0], 2.*b[l], 2.*b[2]) ) 

self . incoming = ( c[l], c[2], c[3], c[0], c[2], c[3], 
c(03, c[l], c[3], c[0], c[l], c(21 ) 

# Construct the coefficients used to form the charged particle 

# in-flow terms, 
self.zz = [] 

self.zz = [ 4.*b[0]', 4.*b[l], 4.*b[2], 4.*b[3] ] 

# 

def rhs(self, geom, qvec, soldata, solved, qcsddata, psiE) : 

"Forms the right hand side of the LD linear tet Boltzmann Equation" 

# Construct the in-flow terms from the adjacent cells. 
XX = [] 

j = -1 

for ivrtx in range (geom, nvrtx) : 

for if ace in range (geom. nf aces) : 
if ( ivrtx == iface ) : 

continue 
else : 

D = j+1 

i = geom. adjvrtx( self . icell, ivrtx, iface) 
if ( i == geom. bdry_f lag) : 

# This demo assumes vacuiun boundary conditions. 

# Other boundary conditions such as reflective, periodic, 

# albedo, or source conditions could be 

# implemented here by setting the value of x 

# appropriately. 
X = 0. 

else: 

# Make sure that no un-set up-wind information 

# is being requested. 

if ( ( not solved[ i[0] ] ) and ( self . incoming [j ] ) ): 
print "cell = self.icell 
print "adj cell= i[0] 

print "Error: attempt to use data from unsolved cell" 
exit < ) 

X = soldata[i[0] ] [i[l] ] 
XX. append (x) 
self.bevc = [] 
for i in range (geom. nvrtx) : 
X = 0. 

for j in rsuige ( 12 ) : 

X = X + self .yy[i] [j]*xx[j] 
self . bvec . append ( x ) 
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# If this is a charged particle problem, create the terms 

# arising from the cyclic reduction of the CSD operator, 
if (self . icp) : 

5 qcsd = qcsddata[self -kk] [self . icell] 

self .gamma = 0. 

for iface in range (geom. nf aces) : 

i = geom. adj cell (self . icell, iface) [0] 
if ( i == geom.bdry_f lag[0] ) : 
10 continue 

else: 

self. gamma = self. gamma + psiE[i] *self . zz [if ace] 
self .sum_qcsd = qcsd [0] +qcsd [1] +qcsd [2 } +qcsd[3 ] 
XX = self .kl*self .k2* (self . sum_qcsd + self. gamma) 
15 cpO = XX + 0.2* (2 , *qcsd[0] +qcsd[l] +qcsd[2]+qcsd[3] ) 

cpl = XX + 0 . 2* (2 . *qcsd[l] +qcsd[0] +qcsd[2] +qcsd[3] ) 
cp2 = XX + 0.2* (2 . *qcsd[2]+qcsd[l] +qcsd[0]+qcsd[3] ) 
cp3 = XX + 0.2* (2.*qcsd[3]+qcsd[l] +qcsd[2]+qcsd[0] ) 
else: 

20 cpO=0.; cpl=0.; cp2 =0.; cp3 = 0. 

# Add in the fixed source, scattering source, CSD terms, et al. 
r = geom. vol (self . icell) /20 . 

self.bvec[0] = self.bvec[0] + r* (2 . *qvec [0] + qvec[l] + 
25 qvec[2] + qvec[3]) + cpO 

self.bvec[l] = self.bvec[l] + r* (2 . *qvec [1] + qvec[0] + 

qvec[2] + qvec[3]) + cpl 
self.bvec[2] = self.bvec[2] + r* (2 . *qvec [2 ] + qvec[l] + 

qvec[0] + qvec[3]) + cp2 
30 self.bvec(3] = self.bvec[3] + r* (2 . *qvec [3 ] + qvec[l] + 

qvec[23 + qvec[0]) + cp3 



35 def ludecomp (self ) : 

"Performs lu decomposition in place using Grout's Method" 

# This method would likely be replaced in a production code by 

# a call to a highly optimized linear algebra library 
40 # function. 

itmp = [] 
self.ib = [] 
for i in range (self .n) : 
45 itmp.append(i) 

self . ib . append ( i ) 



50 



# Loop over columns 

for j in range ( self. n-1) 



# Find pivot 
imax = j 

colmx = abs ( self . amat [ j ] [ j ] ) 
for i in range (j +1 , self . n) : 
55 dum = abs ( self . amat [i] [j ] ) 

if (dum > colmx) : 
imax = i 
colmx = diim 

60 # Pivot this row, if required 

i f ( imeuc ! = j ) : 

for jj in reLnge(self .n) : 

dum = self .amat [imax] [jj ] 
self .amat [imax] [jj ] = self .amat [j ] [jj ] 
65 self .amat [j ] [jj ] = dum 

dum = itmp [imax] 
itmp [imax] = itmp[j] 
itmp[j] = dum 
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# Perform Crout ' s method algebra 
if { self -amat [j] [j] ==0. ): 

print "Error: zero pivot encountered in ludecomp" 

exit () 

self .amat [j ] [j ] = 1 . /self . amat [j ] [j ] 
for i in range ( j +1 , self . n) : 

self . amat [i] tj ] = self . sonat [i] [j ] *self .amat [j ] [j ] 
for jj in range (j+1) : 

for ii in range (j j+1, self .n) : 

self .amat [ii] [j+1] = self . amat [ii] [j+1 1 - \ 

self .amat [ii] [j j] *self .amat [ j j ] [j+1] 

# Invert the last diagonal element 
ii = self.n-1 

if (self .amat [ii] [ii] == 0.): 

print "Error: zero pivot encountered in ludecoit^" 
exit ( ) 

self . amat [ii] [ii] = 1 . /self . amat [ii] [ii] 



# Store the permutation vector 
for i in range ( self . n) : 
self .ib[itmp[i] ] = i 

# 



def fbsub(self ) : 

"Given a lu decomp cuid a source, solves the system via fwd cuid back 
subs t i tu t ion " 



# This method would lilcely be replaced in a production code by 

# a call to a highly optimized linear algebra library 

# function. 



# permute the source 
X = self.bvec[:] 

for i in range ( self . n) : 

self.bvec[ self.ib[i] ] = x[i] 

# Forward substitution 
for j in range (self . n-1) : 

for i in rajige (j+1, self .n) : 

self.bvec[i] = self.bvec[i] - self . amat [i] [j ] *self .bvec [j ] 

# BacJc substitution 

for j in range (self .n-1, -1, -1) : 

self.bvec[j] = self .bvec[j] *self .amat [j] [j] 
for i in range( j-1, -1, -1) : 

self .bvec [i] = self . bvec [i] - \ 

self .amat [i] [ j ] *self .bvec [ j ] 

# If this is a charged particle problem, then compute 

# the energy slope term in the solution and create the 

# CSD source term for the next group, 
if (self . icp) : 

self.psiE = self .k2* (self . sum_qcsd + self .gamma - self.kl*( 

self .bvec [0] + self .bvec [1] + self. bvec [2] + self .bvec [3] ) ) 

self.qcsd = [ self . kl* (self .bvec [0] - self.psiE), 
self .kl* (self .bvec[l] - self.psiE), 
self .kl* (self .bvec[2] - self.psiE), 
self .kl* (self .bvec[3] - self.psiE) ] 



# 

# end of bte.py 

# 

# 

# 

# File: edit.py 
# 

# Defines procedures for producing edits. 
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# 

# Once a solution is calculated, links to any number of various 

# post-processing utilities can be generated easily. A sample 

# visualization link file generator is included here for 

# debug €uid demonstration purposes. 
# 

# Author: John McGhee, Radion Technologies 
# 

# Date : 19 Feb 2004 

# 

# 

# $Id: edit.py,v 1.6 2004/03/05 21:51:31 mcghee Exp $ 
# 

# 

def gmvlink (geom, momdata, outfile) : 

"Writes GMV format visualization link file" 

# GMV is a publicly available visualization utility 

# that can be obtained from Los Alamos National Laboratory, 

# Los Alamos, NM. http://www-xdiv.lanl.gov/XCM/gmv/ 

from string import join 
from time import ctime, time 
from math import sqrt 

# Write the GMV visualization file to disk. 

print "\nWriting GMV format visualization link file to disk 

# Open the output file, 
f = open(outfile, "w" ) 

# Write magic cookie. 

f. write ( "gmvinput ascii \n") 

# Write the nodes block 
data__per_line = 10 

f .write (" \nnodes " + str (geom.nnodes) + "\n") 
for j in range (geom. ndim) : 
ic = 0 

for i in geom . node_coord : 
ic = ic+1 

f .write("%e " % i[j] ) 
if (ic % data_per_line == 0) : 
f .write("\n") 
if (ic % data_per_line != 0): 
f .write ("Nn" ) 

# Write the cells block 

f .write (" \ncells " + str (geom. ncells) + "\n" ) 
for i in geom. eel l_,node: 
record = " " 

f ,write( "tet 4\n" ) 

j =[] 

for k in i : 

j . append ( k+1 ) 
record = join(map(str , j ) ) + "\n" 
f . write ( record) 

# Write the material block 

f .write ( "\nmaterial" + " " + str (geom. nregions) + " 0\n") 
tmp=[] 

for j in range (geom. nregions ) : 

tmp. append ( "mat #- " + str(j)) 
f .write (join (tmp) + "\n") 
ic = 0 

for j in geom. eel l_region: 
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ic = ic+1 
k = j+1 

f.write("%i " % k) 
if (ic % data_per_line == 0) : 
5 f .write (" \n" ) 

if (ic % data_per_line != 0) : 
f .write("\n'') 

# Write the variable block. Variables printed here are fluxes cuid 
10 # currents for all energy groups. Other variables of specific 

# interest to a particular problem or user can easily be included if 

# required. 

moms = (0, 2, 3, 1) # reorder currents to x,y,z 
f .write (" NnvariableXn" ) 
1 5 nmom = min ( 4 , len (momdata) ) 

ngroups = len (momdata [ 0 ] ) 
for ig in range (ngroups) : 
if (ngroups > 1) : 

chdum="-g" + str(ig) 
20 else: 

chdum= " " 
for ii in range (nmom) : 
imom = moms[ii] 
if (imom == 0) : 
25 xnorm = 1, 

f .write (''phiO%s 0 \n" % chdum) 
elif (imom == 1) : 

xnorm = l./sqrt(3.) 
f .write ( "Jz%s 0 \n" % chdum) 
30 elif (imom ==2): 

xnorm = l./sqrt(3.) 
f .write ( "Jx%s 0 \n" % chdum) 
else: 

xnorm = l./sqrt(3.) 
35 f. write ( "Jy%s 0 \n" % chdum) 

ic = 0 

for i in range (geom.ncells) : 
ic = ic+1 
XX = 0 . 

40 j = momdata [imom] [ig] [i] 

kk = 0 
for k in j : 

kk = kk+1 
XX = XX + k 

45 XX = xnorm*xx/ float (kk) 

f.write("%e " % xx) 
if {ic % data_per_line == 0) : 
f .write{ " Nn" ) 
if (ic % data_per_line != 0): 
50 f .write( "\n" ) 

f .write ( " \nendvars\n" ) 



55 



# Write the end of file keyword 
f .write (" \nendgrav NnXn") 

# Close the files, 
f . close ( ) 



60 # Echo success . 

print "GMV visualization link written to file \''%s\''" % outfile 
t = ctime(time( ) ) 
print t, "\n" 

65 

# 

# end of edit.py 

# 
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Scmiple of standard out from the Frost transport solver. 
$Id: frost. out_std,v 1.4 2004/03/10 04:24:22 mcghee Exp $ 



FROST - FiRst Order Sn Transport 
Radion Technologies 

Mesh data successfully read from file mesh.inp 
ncells = 862 

nnodes = 207 

nbdry_faces= 200 

Calculation of mesh geometry data complete, checking results.... 
Mesh data checks passed. 

Problem Description 
ngroups= 2 
nang = 18 

nmat = 2 
nscxs = 2 
nmom = 4 

Neutral Particle Problem 
Calculating solution in group 0 



0 


delta= 


420.854728454 




1 


delta= 


3. 


.38884624607 




2 


delta= 


0. 


.027504091683 




3 


delta= 


0. 


.000224887198983 


4 


delta= 


1. 


.85051804069e- 


06 


5 


delta= 


1. 


.530461755196- 


08 


6 


delta= 


1. 


.270646961586- 


10 


7 


delta= 


1. 


.057961226936- 


12 


8 


delta= 


8. 


.827322943176- 


15 



Group 0 converged . 



Calculating solution in group 1 

0 delta= 0.40013959077 

1 delta= 0.00127508909827 

2 delta= 4 . 07191272391e-06 

3 delta= 1. 30328388021e-08 

4 delta= 4 . 18163292129e-ll 

5 delta= 1 . 3448187674e-13 
Group 1 converged. 



Writing GMV format visualization link file to disk 

GMV visualization link written to file "gmv.out" 
Fri Mar 5 14:39:04 2004 

Dose for entire mesh: 0.00540526704499 

Results in group 0, moment 0, (x=0.5, y=0.6, z=0.7) 
phi0= 0.49361140853 , in cell number: 646 

Balance Table — 

leakage = [0.71825356732193335, 0.018107227200817064] 

other removal = [0.2817464257358761, 0.019458961361966235] 

total source = [0.99999999999999756, 0.037566190098116937] 

balance = [6.9421881576658961e-09, 4 , 0870091755351723e-08] 



Transport Solution Complete. 
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# ! /usr /bin /python 
# 

# 

# File: frost. py 



# 

# FROST - FiRst Order Sn Transport . 
# 

# A simple example of a 3D numerical solver for the first order 

# form of the Boltzmann Transport 

# Equation (BTE) using a linear discontinuous finite element 

# spatial discretization, a discrete ordinates angular 

# discretization, and a multi-group energy discretization. 

# Written in Python { http: //http: //www. python. org/ ). 

# Anisotropic scattering, multiple-material, multiple-particle, 

# and charged-particle problem types are supported. Also 

# provides an option for a first scattered distributed source 

# solution demo mode. Any type of neutral or charged particle may be 

# transported including neutrons, photons, protons, and electrons. 
# 

# This implementation is not optimized for computational 

# efficiency but rather seeks to provide a relatively simple and 

# clear example of the principles composing the algorithm. This 

# demo is based on a tetrahedral finite element mesh for 

# purposes of simplification, but the 

# techniques used herein can and have been applied on other mesh 

# types as well. Like-wise the solution is assiuned to be 

# linear within a tet, but higher order solution trial spaces 

# can be applied as well in a stemdard finite element manner. 
# 

# Method of Solution: 

# The method of solution demonstrated herein consists of 

# ordering the cell equations into 

# a block lower triangular form, sweeping the mesh for each 

# discrete ordinate, and then computing the scattering source via 

# source iteration. Other solution methods are possible. Of note 

# is recent work at Los Alamos on Krylov based solvers rather 

# than source iteration. Various parallelization opportunities 

# are associated with these solvers that may make them an 

# attractive alternative in the future. 
# 

# Source code size and performance: 

# This demo code is approximately 2,400 lines. It runs the 

# included demo problem in approximately 160 seconds on a 2.5 

# Mhz desktop PC. A Fortran code in production use which 

# provides more efficiency, flexibility, error checks, a user 

# friendly interface, etc. is approximately 

# 50,000 lines and performs the same calculation in a fraction of 

# second . 
# 

# Parallelization: 

# This is a serial code. Various methods for parallelization are known 

# and have been implemented in production versions of this 

# algorithm. In general parallelization does not substantially 

# alter the structure of the algorithm and we elect not to 

# include parallelization herein in the interest of clarity. 

# However, comments are inserted where appropriate to indicate 

# where opportunities for parallelization exist. 
# 

# Disclaimer: 

# The purpose of this code is to demonstrate the essential 

# characteristics of the solution algorithm. It is not a 

# production ready solver. As such, the code has been subjected 

# to a minimal set of testing and debug procedures. The results 

# compare favorably with other more well tested solvers. It produces 

# qualitatively reasonable results and passes fundamental tests 

# for accuracy such as particle balance. However, it has not 
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# been tested extensively or thoroughly debugged on a wide 

# variety of problems. Caveat Emptor. 
# 

# A note on Python syntax: 

5 # Indentation is significant ! Indentation is used to denote 

# different blocks of code such as bodies of functions, 

# conditionals, loops, and classes. Care should be exercised 

# when copying or reformatting Python source code to avoid 

# unintentionally altering the functionality of 
10 # the code. 

# 

# Usage: ./frost.py [fsds] 
# 

# The optional argument "fsds" triggers a first-scattered 
15 # distributed source solution mode option. 

# 

# Input is read in from three files: 

# — aquad.inp, angular quadrature data. 

# — mesh.inp, problem geometry and computational mesh. 
20 # matprop.inp, material properties 

# 

# Fixed sources are defined in a the file "fsrc.py". Two 

# sources are provided, both isotropic in group 0 with 

# strength of 1. One source is uniformly distributed, and 
25 # the other is a point source located at (0,0,0). 

# 

# Author: John McGhee, Radion Technologies 

# Date : 19 Feb 2004 
# 

30 # . 

# 

# $Id: frost. py,v 1.41 2004/03/10 17:17:26 mcghee Exp $ 
# 

# 

35 

# Import methods and classes used to solve the transport equation, 
from sys import exit, argv 

from geom import Geometry 

from bte import BteEquation 
40 from acjuad import Quadrature 

from fsrc import FixedSource 

from matp import MaterialProps 

from sord import SolutionOrder 

from edit import gmvlink 
45 from scat import self_scatter , in_scatter 

from fsds import Fsds 

from trace import TetTrace 

# Write a header. 

50 print "XnFROST - FiRst Order Sn Transport" 
print " Radion Technologies\n" 

# Get the command line parameters 
i = len(argv) 

55 if (i == 1) : 

fsds_prob = 0 
elif ( i == 2 ) : 

if (argv[l] == "fsds"): 
fsds_prob = 1 
60 else: 

print "Error: unrecognized command line parameter \ " " , argv [1] , " \ " 
print "Valid options are: no-argrument | fsds" 
print 
exit ( ) 
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# Set convergence criteria cuid max allowed iterations for the source 

# iteration process. Ordinarily this would be set by the user on a by 

# problem basis. 
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eps = l.e-12 

iter^max = 20 

# Read mesh geometry. Tets are assumed here, but 

# other geometries are possible. Higher order solution trial 

# spaces (i.e. sub-parametric, or p-ref inement ) can also 

# be used in the standard finite element manner, 
meshdata = Geometry ( "mesh. inp" ) 

# Read angular Quadrature data. Any angular quadrature 

# can be used depending on the requirements of the problem at hand. 

# For example, Radion has developed special 3D lobatto-chebychev sets 

# for use with problems which contain pleuie wave sources, 
qdata =Quadrature ( " aquad . inp " ) 

# Read multi-group material properties. Multiple material problems 

# with anisotropic up, down, and wi thin-group scatter are supported. 

# Scattering cross sections are provided in terms of an expansion in 

# Legendre coefficients. Multiple particle type problems can be handled 

# transparently by the multi -group method by simply assigning 

# different groups to different particle types. For this demo code 

# there is assumed to be a one-to-one 

# correspondence between the materials and the mesh regions, a 

# production version would provide more flexibility in material 

# assignment. Also of note, adjoint solutions to the transport 

# ec[uation are easily provided for by performing simple modifications 

# to the material properties and their indexes. There are normally 

# only a few tens of materials in a transport problem, and many 

# thousands of computational cells. ConsidereJ^le memory savings can 

# be achieved by storing the material properties by material 

# and assigning material properties to cells only on an as -needed 

# basis. 

matprop = MaterialProps ( "matprop . inp" ) 

# Setup fixed volume source. Full anisotropic, spatially varying, multi-group 

# sources are supported. For this demo code the fixed source values 

# are provided by a pre-defined function. In a production code, 

# provisions would be made to allow the user to specify the fixed 

# source characteristics in detail as part of the problem setup, 
fsdata = FixedSource( qdata, meshdata, matprop) 

# Setup fixed boundary source . 

# This demo assumes volume sources only, but boundary sources 

# can be handled in a similar meoiner to volume sources. Botindary 

# sources are entered into the solution algorithm in the rhs method of 

# the BteEquation class. 

# Perform cuialytic un-collided flvix calculation. 

# If desired an un-collided solution component can be calculated using 

# analytic, high Sn order, or other techniques. This un-collided 

# component can then be used to compute a first collided source to 

# initiate a follow on calculation to complete the solution. Among 

# other things this is often useful as a ray-effect mitigation 

# technicjue. Second and/or subsequent collided solution components and 

# their associated scattering sources can also be 

# calculated and combined in a similar manner if desired, 
if (fsds prob) : 

if (matprop. icp) : 

print "Error: fsds option is not supported for charged particle problems" 

exitO 

fsdsdata = Fsds (matprop, meshdata, qdata, fsdata) 

# Create an array for storing the euigular source, 
qvec = [ 0., 0., 0., 0. ] 

# Create angular solution storage array, 
soldata = [] 



Docket No. 200313554-1 



74 



for i in xrange (meshdata .ncells) : 

soldata. append ( [0 - , 0., 0., 0.]) 

# Create the moment solution storage array. For a reduction in memory 
5 # at the expense of runtime, this array (and others) could be stored 

# on disk and read into memory group -by- group as needed. This array 

# can be initialized from a previous solution stored on disk, which 

# can save substantial computational work in the iterative solution 

# process by providing a good starting guess. This is especially 

10 # advantageous in problems such as perturbation analysis, where the 

# solution is changing minimally from one analysis to the next, 
momdata = [] 

for k in range (qdata. nmom) : 
momdata . append ( [ ] ) 
1 5 for j in range (matprop . ngroups ) : 

x= [] 

for i in xrange (meshdata .ncells) : 

X . append {[0., 0., 0., 0.]) 
momdata [k] . append (x) 

20 

# Create a tmp storage array for source iteration 
oldmom = [ ] 

for k in range (qdata . nmom) : 
oldmom . append ( [ ] ) 
25 for j in range (matprop. ngroups) : 

x= [] 

for i in xrange (meshdata. ncells) : 

x. append ( [0. , 0., 0., 0.]) 
oldmom[k] .append(x) 

30 

# If this is a charged particle problem, then the transport equation 

# is altered to include the continuous-slowing-down (CSD) term from 

# the Boltzmann-Fokker-Planck operator. We difference this term using 

# a linear discontinuous treatment in energy which introduces an 

35 # additional unknown into the solution trial space. Additional storage 

# is required for the treatment of this unknown, which we set up below. 
psiE = [] 

qcsd = [ [], [] ] 
csdptr_g = 1 
40 csdptr_gml = 0 

if (matprop. icp) : 

for k in range (2) : 

for j in range (qdata. nang) : 

^= " 

45 for i in xreuige (meshdata. ncells) : 

X. append ( [0., 0., 0., 0.] ) 
qcsd[k] . append (x) 
for i in xrange (meshdata .ncells) : 
psiE . append ( 0 . ) 

50 

# Create storage arrays for a balance table, 
abs = [] 

leakage = [ ] 
source = [ ] 
55 balance = [] 
csd_src = [] 
csd_rem = [ ] 

for j in range (matprop. ngroups ) : 
abs . append ( 0 . ) 
60 leakage . append ( 0 . ) 

source . append ( 0 . } 
balance . append ( 0 . ) 
cs4_src . append ( 0 . ) 
csd_rem . append ( 0 . ) 
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# Determine sweep ordering. For a reduction in memory at a cost of 

# increased runtime, this could be computed on the fly for each angle 

# inside of the source iteration loop. Also, a fine grained 
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# parallelization can be accomplished at this step by assignment of 

# cells to a multiprocessor solution schedule based on a careful analysis 

# of the graph associated with the problem mesh. This technique has been 

# implemented in research and development prototypes, 
swpdata = SolutionOrder (meshdata, qdata) 

# Echo a few details to std out to identify the problem parameters, 
print "Problem Description" 

print " ngroupss " , matprop . ngroups 
print " nang = qdata.nang 
print " nmat = matprop. nmat 
print " nscxs = matprop . nscxs 
print " nmom = qdata.nmom 
if (matprop. icp) : 

print " Charged Particle Problem" 
else: 

print " Neutral Particle Problem" 
if (fsds^^rob) : 

print " First Scattered Distributed Source Problem" 
print " nrays = f sdsdata . nrays 

print " average cells per ray = %5.1f" % f sdsdata. ncells_per_r ay 
print 

# If upscatter is present, then an additional "outer" iteration would 

# be wrapped around the energy group loop at this point. This demo 

# assumes downscatter only, so the outer iteration loop is not 

# present. Also, for some classes of multi-particle problems it may be 

# advantageous to solve the block of groups associated with each 

# particle as a separate calculation in order to optimize the 

# computation for each particle type. 

# Loop over energy groups from high energy to low. We use the nuclear 

# engineering convention for group numbers (ie. group 0 is the highest 

# energy group, and group "ngroups" is the lowest energy group, 
for ig in range (matprop. ngroups) : 

print "Calculating solution in group ", ig 

delta = l.# a measure of the change in the solution from one 

# iteration to the next 
icoTint = -1 # the iteration counter 

# For charged particle problems there is a source from the CSD 

# term from the group above in energy. To save memory, we only 

# store the source from the preceding group and the current 

# group. Here we swap the pointers for the current and 

# preceding groups. 

X = csdptr_gml 

csdptr_gml = csdptr_g 
csdptr_g = x 

# loop over the wi thin-group scattering (self -scatter) source 

# until converged for this group, 
while (delta > eps) : 

# Clear the balance table 
csd_src [ig] = 0 . 
csd_rem[ig] = 0 . 
leakage [ig] = 0. 
abs[ig] =0. 

source [ig] = 0. 

# increment the iteration counter 
i count = i count + 1 

# trap a likely error condition to prevent infinite iteration, 
if (icount > iter_meuc) : 

print "Error: maximum allowed number of inner iterations exceeded. 



print "Iterations 
print "Current delta 
print "Convergence crit 
break 



icount 

delta 

eps 
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# save the old solution for comparison with the result of this 

# source iteration step. 

for imozn in range (qdata. nmom) : 

for i in xrange (meshdata. ncells) : 

for j in xrange (meshdata. nvrtx) : 

oldinom[imom] [ig] [i] [j] = momdata [imom] [ig] [i] [jl 
momdata[imom] [ig] [i] [jl = 0. 

# loop over angles. Angles proceed independently in the 

# absence of implicit boundary conditions. A coarse grained 

# parallelization can be implemented by assigning individual 

# eungles to separate processors at this point, 
for kk in range ( gdata . nang } : 

solved = [] 

for ic in xrange (meshdata. ncells ) : 

solved . append ( 0 ) 
omega = gdata. qdpt (kk) 
swp = swpdata. order (kk) 
m2d = gdata. mom2dis (kk) 



# loop over cells 

for ic in xrange (meshdata . ncells) : 

# find next cell to solve 
icell = sv7p[ic] 

# create BTE for this cell 

xstot = matprop. sigt ( ig, meshdata. region (icell) ) 
beta = matprop . sigstp ( ig, meshdata. region (icell) ) 

# create the source moments, note that the fixed and 

# inscatter sources could be calculated and stored outside the 

# self -scatter loop for increased computational efficiency, 
if (fsds_prob) : 

fixed = fsdsdata. fsdsMom ( icell, ig) 
else: 

fixed = fsdata. fsrc (icell , ig) 
inset = in_scatter (icell, ig, meshdata, matprop, oldmom, qdata) 
self set = self_scatter (icell, ig, meshdata, matprop, oldmom, gdata) 

# convert from moment trial space to discrete trial space 
for j in range (meshdata . nvrtx) : 

qvec[jl = 0. 

for i in range (qdata . nmom) : 

qvec[j] = gvec[j] + m2d[i]*( 

fixed[i][j] + insct[i][j] + self set [i] [j ] ) 

# solve BTE for this cell 

b = BteEguation (meshdata, icell, matprop. icp, omega, 
xstot, beta, kk) 

b . ludecomp ( ) 

# for increased speed at the cost of more memory, the 

# ludecomp could be stored here for use in future iterations, 
b. rhs (meshdata, qvec, soldata, solved, 

qcsd [csdptr_gml] , psiE) 

b. fbsub( ) 

solved [icell] = 1 

# store the solution 
soldata [icell] = b.bvec[:] 
if (matprop. icp) : 

psiE[icell] = b.psiE 

qcsd[csdptr_g] [kk] [icell] = b.qcsd[:] 

# return for next cell in this cuigle. 

# accumulate the moments of the angular solution for 

# computation of scattering source and edits, 
for m in range (gdata. nmom) : 

x = gdata. dis2mom(kk) 

for j in range (meshdata. ncells) : 

for i in range (meshdata. nvrtx) : 

momdata[m] [ig] [j] [i] = momdata [m] [ig] [j] [i] + \ 

soldata[j] [i]*x[m] 
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# accumulate leakage for the balance table edit 
for bf in meshdata.bdry_f ace : 

XX = 0.; icell = bf[0]; iface = bf[l] 
yy = meshdata. avect (icell, iface) 
5 for j in range (meshdata.ndim) : 

XX = XX + oinega[ j ] *yy [ j ] 
if (xx > meshdata.dplimit) : 

xO = meshdata. f ace_node [if ace] [0] 
xl = meshdata. face_node [iface] [1] 
10 , x2 = meshdata, face_node [iface] [2] 

leakage [ig] = leakage [ig] + \ 

( 1 . /3 . ) * omega [meshdata . ndim] *xx* ( \ 
soldata [icell] [xO] + \ 
soldata [icell] [xl] + \ 
15 soldata[icell] [x2] ) 

# accumulate the CSD energy source and removal terms for 

# the balance table, 
if (matprop . icp) : 

20 wgt = omega [meshdata . ndim] ; xO = 0 . ; xl = 0 . 

for j in range (meshdata. ncells) : 

for i in range (meshdata .nvrtx) : 

xO = xO + qcsd[csdptr_gml] [kk] [j] [i] 
xl = xl + <3csd[csdptr_g] [kk] [j] [i] 
25 csd_src[ig] = csd_src[ig] + wgt*xO 

csd_rem[ig] = csd^rem[ig] + wgt*xl 

# return for next angle in this energy group 

30 # Compute a measure of the change in the solution 

delta = 0. 

for i in xrange (meshdata .ncells) : 

for j in xrange (meshdata. nvrtx) : 

delta = delta + (oldmom[0] [ig] [i] [j] - momdata [0] [ig] [i] [ j ] ) **2 
35 print " icount, "delta= delta 

# At this point convergence of the scattering source can be 

# accelerated with any number of pre-conditioning techniques 

# such as diffusion synthetic acceleration (DSA) . This is 

40 # essential for problems with a high within groups scattering 

# ratio. 

# return for next wi thin-group scattering source iteration. 

45 print "Group ig, "converged." 

print 

# Confute other balance table entries 
imom = 0 

50 for j in range (meshdata. ncells ) : 

if (f sds _prob) : 

XX = f sdsdata. f sdsMom( j , ig) [0] 

else: 

XX = f sdata. fsrc( j , ig) [0] 
55 yy = in_scatter ( j , ig, meshdata, matprop, momdata, qdata) [0] 

for i in range (meshdata . nvrtx) : 

abs[ig] = abs[ig] + momdata [imom] [ig] [ j ][ i ] * ( \ 

matprop. sigt ( ig, meshdata, region ( j ) ) - \ 
matprop. sigsct ( 0, ig, ig, meshdata. region(j ) ) )* \ 
60 meshdata . vol ( j ) 

source [ig] = source [ig] + meshdata. vol (j )* ( xx[i] + yy[i] ) 
abs[ig] = 0.25*abs[ig] + csd_rem[ig] 
source [ig] = 0 . 25*source [ig] + csd_src[ig] 

65 # return for next energy group 

# At the end of the loop over energy groups, if up-scatter is present, 

# diffusion synthetic or other acceleration technic[ues similar to 
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# those employed for the acceleration of the self scattering source 

# can be employed to speed convergence at this point. Since this demo 

# is restricted to down-scatter for purposes of simplicity, 

# no outer loop acceleration algorithms are necessary. 

# If this is a first scattered distributed source problem, then we 

# must add in the un- collided solution component to the collided 

# component to complete the solution, 
if (fsds_prob) : 

f sdsdata . total_solution (momdata) 

# Solution is contplete at this point. Now we begin post-processing. 

# If the solution is saved to disk here, then any desired post processing 

# subsequent to this run can be accomplished without the work of 

# re-computing the solution. 

# Write a visualization link file, 
gmvlink (meshdata , momdata , ■ gmv . out " ) 

# Compute an example dose edit by integrating the results over 

# the finite element mesh. Other edits can be computed in a 

# similar manner as required. If edits at an arbitrary point are 

# required this is also easily computed as the finite element trial 

# space provides a rigorous means of interpolation at any point 

# in the mesh. 

for ig in range (matprop.ngroups) : 
result = 0 . 
imom - 0 

for j in range (meshdata. ncel Is) : 

for i in range (meshdata. nvrtx) : 

result = result + momdata [imom] [ig] [j ] [il * \ 

matprop.sigreac ( ig, meshdatja. region( j ) )* \ 
meshdata . vol ( j ) 

result = 0.25*result 

print "Dose for entire mesh: result 

print 

# Compute an example result at an arbitrary point in the mesh. Call 

# the TetTrace interpolation routines to produce the correct LDFEM 

# value for group 0, moment 0 at point (0.5, 0.6, 0.7). Also find the 

# cell number containing this point for reference, 
pt = (0.5, 0.6, 0.7) 

trace = TetTrace (meshdata) 

X = trace . fieldVal ( pt, momdata [0] [0] ) 

i = trace. eel lNumb( pt ) 

print "Results in group 0, moment 0, (x=0.5, y=0.6, z=0.7)" 

print " phi0= x, in cell number: i 

print 

# Complete a balance table and print. The balance table can be used 

# as an error checking tool, a convergence metric, or as 

# another set of edit quantities of interest to the user, 
print "Balance Table--" 

print " leakage = " , leakage 

print " other removal = " , abs 
print " total source = " , source 
for ig in range (matprop . ngroups ) : 

if (source[ig] != 0.): 

balcince[ig] = 1. - (abs [ig] +leakage [ig] ) /source [ig] 

else: 

balcuice[ig3 = abs [ig] +leakage [ig] 
print " balance = " , balance 

# Post processing complete. 

# Write a footer and halt. 

print "XnTransport Solution Complete. NnNn" 
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# 

# end of frost.py 

# 

5 # 

# 

# File: fsds.py 
# 

# Provides information related to the 

10 # "First-Scattered- Distributed- Source" solution technique. 

# 

# Computes the un-collided component of a treuisport solution 

# using an emalytic ray tracing technique, then forms a first 

# scattered distributed source. 

15 # 

# This siznple exan^le assumes a single isotropic point source. 

# A general purpose routine would provide for 

# the entry of any number of user specified multi-group 

# anisotropic sources . 

20 # 

# Author: John McGhee, Radion Technologies 
# 

# Date : 19 Feb 2004 

25 # 

# 



30 



35 



50 



55 



# $Id: fsds.py, V 1.9 2004/03/10 17:17:26 mcghee Exp $ 
# 

# 

from sys inqport exit 
class Fsds : 

"First Scattered Distributed Source methods and data" 



def init (self, matprop, meshdata, qdata, fsdata) : 

"Computes uncollided solution component and first scattered source" 
from math import sqrt, pi, exp 
from sph import sphfun 
40 from trace import TetTrace 

self.ngroups = matprop. ngroups 

self.nvrtx = meshdata . nvrtx 

self.nmom = qdata. nmom 

45 self .ncells = meshdata .ncel Is 

self .ncells_per_ray = 0 

self.nrays = 0 

self .uncolMom = [] #[icell] [ig] [imom] [ivrtx] 

self.sctsrc = [] # [icell] [ig] [imom] [ivrtx] 



# A production algorithm would likely provide for multiple 

# sources by including a loop over all the sources in the 

# problem at this point. As we assijme only one source for this 

# demo, no loop over sources is required. 



# Get the point source location and strength. An isotropic 

# source is assumed herein for simplicity, but in general this 

# is not a restriction. In a production code, the source 

# details would be specified by the user at problem setup. 
60 srcID = 0 

(ray_starting_point, ptSource) = fsdata. isopsrc ( srcID) 

# To form the uncollided solution, we track 

# from the source to the Gaussian integration point (s) on each 
65 # cell. The first scattered source is then rigorously computed 

# by using finite element integration rules on the cell. Here 

# we go. 
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# As an example we use a four point quadratic Gaussian 

# integration rule on linear tetrahedra. Other integration 

# rules could be used if desired, 
self.nqdpt = 4 

5 alpha = 0.05* (5. + 3.*sqrt(5.) ) 

beta = 0.05* (5. - sqrt{5.) ) 
gamma = 4 , *alpha - 3 . *beta 
delta = 4.*beta - {alpha + 2.*beta) 

10 # Loop over cells computing the un-collided solution component and 

# the first scattered source, 
ncpr = 0 

for icell in xrange(self .ncells) : 

nO = meshdata. eel l_node [icell] [0] 
15 nl = meshdata.cell_node[icell] [1] 

n2 = meshdata. cell_node [icell] [2] 
n3 = meshdata. eel l_node [icell] [3] 
cO = meshdata.node_coord[nO] 
cl = meshdata.node_coord [nl] 
20 c2 = meshdata. node_coord[n2] 

c3 = meshdata .node_coord [n3 ] 
phi= [] #[iqdpt] [ig] [imom] 
# Loop over quadrature points 
for iqdpt in range (self .nqdpt) : 
25 # Compute the coordinates of the Gaussian integration point 

if (iqdpt == 0) : 

ray_ending_point = ( 

alpha*c0[0] + beta* (cl [0] +c2 [0] +c3 [0] ) , 
alpha*c0[l] + beta* (cl [1] +c2 [1] +c3 [1] ) , 
30 alpha*c0[2] + beta* (cl [2] +c2 [2] +c3 [2] ) ) 

elif (iqdpt == 1) : 

ray_ending_point = ( 

alpha*cl[0] + beta* (cO [0] +c2 [0] +c3 [0] ) , 
alpha*cl[l] + beta* (c0[l]+c2 [l]+c3 [1] ) , 
35 alpha*cl[2] + beta* (cO [2] +c2 [2] +c3 [2] ) ) 

elif (iqdpt == 2) : 

ray_ending_point = ( 

alpha*c2[0] + beta* (cO [0] +cl [0] +c3 [0] ) , 
alpha*c2[l] + beta* (cO [1] +cl [1] +c3 [1] ) , 
40 alpha*c2[2] + beta* (cO [2 ] +cl [2] +c3 [2 ] ) ) 

else: 

ray_ending_point = ( 

alpha*c3[0] + beta* (cO [0] +cl [0] +c2 [0] ) , 
alplia*c3[l] + beta*(c0[l]+cl[l]+c2[l] ) , 
45 alpha*c3[2] + beta* {cO [2] +cl [2] +c2 [2] ) ) 

# Compute the direction of particle travel. 

omega = ( ray__ending_point [0] - ray_starting_point [0] , 
ray_ending_point ( 1 ] - ray_starting_jE>oint [1] , 
50 ray_ending__point [2] - ray_starting_point [2] ) 

# At this point the algoritlnn trac]cs through 

# the mesh tet by tet from the starting point to the 

# ending point of the ray, accumulating the decay of the source 

# based on the material assigned to each tet and the path 
55 # length in that tet. Alternatively, the ray could be 

# traced through the elements of any other problem 

# related geometry deemed appropriate, for 

# example a voxel based CT scan. The only output 

# required from the traclcing algorithm is the optical 
60 # path length from source to quadrature point. Methods 

# provided in trace. py (TetTrace class) provide an 

# example of how this can be accomplished on tet meshes, 
trace = TetTrace (meshdata) 

ncpr = ncpr + trace . traclcData (icell, ray_ending_point , 
65 ray_starting_jpoint ) 

# Compute the path length and normalization consteuit 
path = sqrt (omega [0] **2 + omega [1]**2 + omega[2]**2) 
xnorm = 1 . / (4 . *pi*path*path) 
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# Loop over energy groups 
yy = [] #[ngroups] [nmoiti] 

for ig in range (ma tprop.ngroups) : 

# compute the optical path 
optical path = 0. 

for i in trace. path: 

iraat = meshdata . cell_region [i [0] ] 

optical_path = optical_path + i [1] *matprop . sigt (ig, imat ) 

# compute the uncollided angular flux at this point 
psi = xnorm*ptSource[ig] *exp(-optical_path) 

# coic^ute the moments of the uncollided fl\ix 

# at this quadrature point. 
XX = [ ] # [nmom] 

for imom in range (qdat a. nmom) : 
1 = qdata , index_lm [ imom] [ 0 ] 
m = qdata . index_lm [ imom] [ 1 ] 

xx.append( psi*sphfun ( 1, m, omega[0], omega [1] , omega[2]) ) 
yy . append ( xx ) 
phi . append ( yy ) #[nqdpt] [ngroups] [nmom] 

# return for the next quadrature point in this cell 

# Now that we have the uncollided solution moments at all 

# the quadrature points in this cell, we map the solution 

# to the cell vertexes . This is done via a rigorous finite 

# element based formula that preserves as many spatial 

# moments of the solution as possible. This 

# amounts to assigning a specific weighted sunx of the solution 

# at the quadrature points to the cell vertexes. 
zz = [] # [ngroups] [nmom] [nvrtx] 

for ig in range (ma tpr op . ngroups ) : 
yy = [] #[nmom] [nvrtx] 
for imom in range (qdata . nmom) : 
XX = n # [nvrtx] 

XX. append ( gamma*phi[0] [ig] [imom] + delta* ( 

phi[l] [ig] [imom] + phi [2] [ig] [imom] + phi [3] [ig] [imom])) 
XX. append ( gamma*phi[l] [ig] [imom] + delta* ( 

phi [0] [ig] [imom] + phi [2] [ig] [imom] + phi [ 3 ] [ ig] [ imom] )) 
XX, append ( gamma*phi[2] [ig] [imom] + delta* ( 

phi [1] [ig] [imom] + phi [ 0 ][ ig] [ imom] + phi [3 ] [ ig] [ imom] )) 
xx.append( gamma *phi [3 ] [ig] [imom] + delta* ( 

phi[l] [ig] [imom] + phi [2] [ig] [imom] + phi[0] [ig] [imom] ) ) 
yy . append ( XX ) 
zz . append ( yy ) 

self .uncolMom. append (zz) #[ncells] [ngroups] [rnnom] [nvrtx] 

# Now that we have the moments of the uncollided solution 

# component on this cell, we can compute the moments of the 

# first scattered source ( par tides /time -length'^ 3 ) . 
imat = meshdata. region (icell) 

zz = [] # [ngroups] [nmom] [nvrtx] 
for isinlc in range (ma tpr op. ngroups) : 
yy = [] #[nmom] [nvrtx] 
for imom in range (qdata . nmom) : 
1 = qdata . index_lm [ imom] [ 0 ] 
XX = [] # [nvrtx] 

for ivrtx in range (meshdata . nvrtx) : 
x = 0. 

for isrc in range (ma tpr op . ngroups) : 

X = X + self .uncolMom [icell ] [isrc] [imom] [ivrtx] * \ 
matprop. sigsct (1, isrc, isin]c, imat) 
XX. append (x) 
yy. append ( XX ) 
zz . append (yy) 

self .sctsrc. append (zz) #[ncells] [ngroups] [nmom] [nvrtx] 
self.nrays = self .nqdpt*self . ncells 

self .ncells_per_ray = f loat (ncpr) /float (self .nrays) 
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# 

def f sdsMom(self , icell, igroup) : 

"Returns the first scattered distributed source for cell icell, and group 

igroup" 

if (icell < 0) or (icell >= self .ncells) : 

print "Error: icell argument out of range in fsds" 
exit ( ) 

if (igroup < 0) or (igroup >= self .ngroups) : 

print "Error: igroup argument out of remge in fsds" 
exit ( ) 

return self . sctsrc [icell] [igroup] # [nmom] [nvrtx] 

# 

def unco lMom(s el f , icell, igroup) : 

"Returns the moments of the un- collided component of the solution. " 
if (icell < 0) or (icell >= self .ncells) : 

print "Error; icell argument out of range in fsds" 

exit ( ) 

if (igroup < 0) or (igroup >= self .ngroups ) : 

print "Error: igroup argument out of range in fsds" 
exit ( ) 

return self . uncolMom[ icell ] [igroup] # [nmom] [nvrtx] 

# 

def total_solution(self , collided) : 

"Combines the uncollided and collided solution components to form the total 
solution. " 

# collided [nmom] [ngroups] [ncells] [nvrtx] 
for imom in range ( self .nmom) : 

for ig in range (self . ngroups) : 

for icell in xrange (self .ncells) : 

for ivrtx in range (self .nvrtx) : 

collided [ imom] [ig] [icell] [ivrtx] = \ 
self . uncolMom [icell] [ig] [imom] [ivrtx] + \ 
collided [imom] [ig] [icell] [ivrtx] 



# 

# end of fsds.py 

# 

# 

# 

# File: fsrc.py 



# 

# Provides information on a fixed source. This source is defined 

# as part of the problem setup. For simplification, this demo 

# provides a single fixed source definition as a Python function 

# call. 
# 

# A production code would ordinarily provide a means for the 

# user to define the details of the source (s) as part of the problem 

# setup. This could be done in a variety of ways such as reading 

# the source definition from disk or selecting the source 

# definition from a library of pre-defined functions . 
# 



# Author: John McGhee, Radion Technologies 
# 

# Date : 19 Feb 2004 

# 

# 

# $Id: fsrc.py, V 1.5 2004/03/05 21:51:31 mcghee Exp $ 
# 
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# 



£roin sys import exit 



class FixedSource: 

"Provides information on the fixed source." 



def init (self, qdata, meshdata, matprop ) : 

"Initializes the fixed source class" 



self.nmoin = 
self . index_lm = 
self.nvrtx = 
self.ncells = 
self.ngroups = 
self. one = [] 
self. zero = [] 



qdata . nmom 
qdata . indexL_lin 
meshdata . nvrtx 
meshdata . ncells 
matprop . ngroups 



for i in range (self .nvrtx) : 
self . one . append ( 1 . ) 
self . zero . append ( 0 . ) 



self .ptSrcLoc = (0., 0., 0.) 
self .ptSrcVal = [1. , ] 
for i in range (1, self .ngroups) : 
self . ptSrcVal . append ( 0 . ) 



def fsrc(self, icell, igroup) : 

"Returns the fixed source for cell icell, and group igroup" 

# This method returns a uniformly distributed, isotropic 

# source with a strength of 1.0 in group 0. 

if (icell < 0) or (icell >= self.ncells): 

print "Error: icell argument out of range in fsrc" 

exit ( ) 

if (igroup < 0) or (igroup >= self . ngroups ) : 

print "Error: igroup argument out of rsmge in fsrc" 
exit ( ) 

X = [] 

for imom in range (self . nmom) : 
1 = self . index_lm [ imom] [ 0 ] 
if ( 1 == 0 ) cund ( igroup == 0 ) : 

X . append ( sel f . one ) 
else: 

X . append ( self . zero } 
return x # [nmom] [nvrtx] 

def isopsrc (self , srcID) : 

"Returns the location and strength of an isotropic point source" 

# This method returns a isotropic source located at (0,0,0) 

# with a strength of 1.0 in group 0. 



return ( self .ptSrcLoc , self .ptSrcVal ) 



# 

# end of fsrc.py 

# 

# 

# 

# File: geom.py 
# 

# Mesh geometry and related functions. 
# 

# Reads in mesh data from disk: ncells, nvertexes, nregions. 
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# vertex coordinates, cell vertexes, 

# region numbers, boundary conditions, cell types. 

# For simplicity, this demo assumes a 3D linear tetrahdedral 
5 # mesh, but other geometries can be treated in a similar 

# manner. The finite element geometric coefficients provided 

# herein are normally computed by a 

# spatial quadrature integration rule. For linear 

# tetrahedra, it so happens that the results of these 

10 # integrations are conviniently expressed in terms of the 

# tet volumes, face areas, and face normal vectors. Other 

# geometries would commonly provide this data as a matrix of 

# coefficients particular to each cell. 

15 # As a matter of computational efficiency, these "geoemtric" 

# coefficients are pre-computed and stored for each cell. This 

# data is then later combined with the Sn angrular quadrature data 

# and material properties data to compute coefficients specific 

# to each space* angle -energy phase space cell. 



20 



25 



40 



55 



# Finite element geomtry data computed for tets: 

# - cell volumes 

# - cell face areas 

# - outward directed cell face unit normal vectors 



# We also compute pointers to adjacent cells and shared vertexes 

# to complete the geometry description. 
# 

# Author: John McGhee, Radion Technologies 

30 # 

# Date : 19 Feb 2004 

# 

# 

35 # $Id: geom,py,v 1.18 2004/03/05 21:51:31 mcghee Ejqp $ 
# 

# 



from sys import exit 

class Geometry: 

"Provides information on the tetrahedral mesh geometry." 



element_geometry = "tetrahedra" 
45 null_flag = (-999,-999, -999) 

bdry_flag = (-998,-998, -998) 
nvrtx = 4 
nfaces = 4 
ndim = 3 

50 face_node = ( (1,2,3), (0,3,2), (0,1,3), (0,2,1) ) 



etc. 



def init (self ,mesh_f ile) : 

"Reads in the tet mesh geometry from disk and computes volumes, areas, 

from string import split, atof, atoi, strip 
from math import sqrt 



# Define a constant for test 
60 self.dplimit = l.e-15 

# Setup for read of geometry data from disk, 
self . name=mesh_f ile 

f = open (mesh_f ile) 
65 ^ record = f.readline() 



# Read in the node coordinates from disk, 
while (record and strip(record) != "nodes" ): 
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record = f.readline{) 
if (record == " " ) : 

print "Error: unable to find nodes keyword on mesh file %s\n" % 

inesh_f ile 

exit ( ) 
self .node_coord = [] 
record = f . readline ( ) 

while (record and strip (record) != " end_nodes " ) : 
j = split (record) 

self . node^coord . append (map ( atof , j [ 1 : 4 ] ) ) 
record = f. readline () 
if (record == ""): 

print "Error: unable to find end_nodes keyword on mesh file %s\n" % 

mesh_f ile 

exit ( ) 

self.nnodes = len(self .node_coord) 

# Read in the cell nodes and region flag from disk. 

# Note the tet node numbering scheme assumes that local node numbers 

# increase in a clockwise direction when viewed from node 0. 

# Also note that the global node numbering scheme is assumed 

# to be zero based, ie the first node in the node list is node 

# "0". Mesh regions are assumed to be id'd with consecutive 

# integers beginning with 0, 
record = f.readlineO 

while (record and strip (record) != "cells"): 

record = f . readline ( ) 
if (record == ""): 

print "Error: unable to find cells keyword on mesh file %s\n" % 

mesh_f ile 

exit ( ) 
self . eel l_node = [] 
self . cell_region = [] 
record = f.readline() 

while (record and strip (record) != "end_cells" ) : 
j = split (record) 
if (atoi( j [1] ) != 4) : 

print "Error: Non-tetrahedral cell id %s. found on mesh file %s\n" % 
( j [0] , mesh_file) 

exit ( ) 

self .eel l_node. append (map (atoi, j [2 : self . nvrtx+2 ] ) ) 
self . cell_region . append ( atoi ( j [ self . nvrtx+2 ] ) ) 
record = f . readline ( ) 
if (record == ""): 

print "Error: unable to find end^cells keyword on mesh file %s\n" % 

mesh_file 

exit () 

self.ncells = len (self . cell_node) 

self .nregions = max (self . cell_region) +1 

# Report successful read of mesh geometry from disk, 
f . close ( ) 

print "Mesh data successfully read from file %s" % mesh_file 
print " ncells = self.ncells 

print " nnodes = ", self.nnodes 

# Compute outward directed cell face area vectors by vector cross 

# product of cell edge vectors. Compute cell voliune by triple scalar 

# product. 

self .eel l_volume = [] 

self .area_vec tor = [] 

self . face_area = [] 

for nlist in self . cell_node: 

nO = nlist [0] 

nl = nlist [1] 

n2 = nlistf2] 

n3 = nlist [3] 

xO = self .node_coord [nO] [0] 
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yo 




self 
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. node. 


.coord [n3] [1] 


z3 




self. 


. node. 


.coord [n3] [2] 


# face 0 






XX 




( 







10 



< 0.5*( (y2-yl)*(23-zl) - (z2-2l) * (y3-yl) ) ), 
15 ( -0.5* ( (x2-xl)*(z3-zl) - (z2-zl)*(x3-xl) ) ), 

( 0.5*( (x2-xl)*(y3-yl) - (y2-yl) * (x3-xl> ) ) ) 
self . area_vector . append (xx) 

self . face_area. append ( sqrt (xx [0] *xx [ 0] + xx[l]*xxtl] + 

xx[2]*xx[23) ) 

20 # face 1 

XX = ( 

( 0.5*( (y3-y0) * (z2-z0) - { z3-z0 ) * (y2-y0 ) ) ), 
( -0.5*( (x3-x0)*(z2-z0) - (23-zO) * (x2-x0) ) ), 
( 0.5*( (x3--x0) * (y2-y0) - (y3-y0) * (x2-x0) ) ) ) 
25 self .area_vec tor. append (XX ) 

self . face_area. append ( sqrt (xx [0] *xx [0] + xx[l]*xx[l] + 

xx[2]*xx[2J) ) 

# face 2 
XX = ( 

30 ( 0.5*( (yl-y0)*(z3-z0) - (zl-zO) * (y3-y0) ) ), 

( -0.5*( (xl-x0)*(z3-z0) - (zl-zO)*(x3-xO) ) ), 
( 0.5*( (xl-x0)*(y3-y0) - (yl-yO ) * {x3-x0) ) ) ) 
sel f . area_vector . append ( xx ) 

self . face_area. append ( sqrt (xx[0] *xx[0] + xx[l]*xx[l] + 
35 xx[2]*xx[2]) ) 

# face 3 
XX = ( 

( 0.5*( (y2-y0) * (zl-zO) - ( z2-z0 ) * (yl-yO ) ) ), 
( ~0.5*( (x2-x0) * (zl^zO) - (z2-z0) * (xl-xO) ) ), 
40 ( 0.5*( (x2-x0) * (yl-yO) - (y2-y0 ) * (xl-xO) ) ) ) 

self , area_vector . append (xx) 

self . f ace_area . append ( sqrt ( xx[0]*xx[03 + xx[l]*xx[l] + 

xx[2]*xx[2]) ) 

# cell volxiine 

45 self .cell_volume, append ( (-l./3.)*( (x3-x0) *xx[0] + 

(y3-y0)*xx[l] + 
<z3-z0)*xx[2] ) ) 

# For cell i, vertex face k, create a pointer to the adjacent 
50 # cell ia and vertex ja. We construct a list of bins so we don't 

# have to search the whole mesh. Each bin contains a list of 

# the cells that share a particular node, 
self . adjacent = [] 

bin = [] 

j3 for i in xrange ( self.nnodes ): 

bin . append ( [ ] ) 
i = -1 

for nlist in self . eel l_node: 
i = i+1 

60 j = -1 

for k in nlist: 
j = j+1 

bin [k] .append ( (i,j) ) 
for icell in xrange < self .ncel Is) : 
65 for ivrtx in xremge (self . nvrtx) : 

for iface in xrange (self .nf aces) : 
if ( iface == ivrtx ) : 
XX = self .null_f lag 
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else : 

found = 0 

node = ( self .cell_node[icell] [self . face_node [if ace] [0] ] , 
self .cell_node[icell] [self . f ace_node [if ace] [1]], 
self .cell_node[icell] [self . f ace_node [if ace] [2]] ) 
for i in bin[self .cell_node(icell] [ivrtx] ] : 
if (i[0] == icell) : 

continue 
j = self . eel l_node [i [0] ] 
if ( (node[0] in j) and (node[l] in j) 
and (node [2] in j) ): 
if ( not ( j [0] in node ) ) : 

jface = 0 
elif ( not ( in node ) ): 

jface = 1 
elif ( not ( j[2] in node ) ): 

jface = 2 
else: 

jface = 3 
found = 1 

XX = (i[0], i[l], jface) 
break 
if (found != 1) : 

XX = self .bdry_f lag 
self .ad j acent . append (xx) 

# Make a list of the boundary faces 
self .bdry_f ace = [] 

for icell in xran9e(self .ncells) : 

for iface in xrange (self .nf aces) : 

if (self .adjcell (icell, iface) [0] == self .bdry_f lag [ 0] ) : 
self .bdry_f ace. append ( (icell, iface) ) 
self .nbdry_f aces = len ( self .bdry_f ace) 
print " nbdry_faces= " , self .nbdry_f aces 

# Check the calculations, and report success (or failure!). 

print " Calculation of mesh geometry data complete, checking results...." 
if ( self . invariant ( ) ): 

print " Mesh data checks passed. \n" 
else: 

print " Mesh data checks failed !\n" 
exit ( ) 



def vol(self,i): 

"Returns the volxraie of cell i" 
if (i < 0) or (i >= self .ncells) : 

print "Error: argument out of range in geometry .vol" 

exit ( ) 

return self .eel l_volume[i] 



def area ( self , i , j ) : 

"Returns the area in cell i of face j " 

from sys import exit 

if (i < 0) or (i >= self .ncells) : 

print "Error: argument i out of range in geometry . area" 

exit () 

if (j < 0) or (j >= self .nfaces) : 

print "Error: argument j out of range in geometry . area " 

exit () 
ii = i*self .nfaces + j 
return self . f ace_area [ ii ] 
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def region(self , i ) : 

"Returns the mesh region number of cell i" 
if (i < 0) or (i >= self .ncells) ; 

print "Error: argument out of range in geometry . region" 
exit () 

return self . cellar egion [i] 

# 

def avect (self / i, j ) : 

"Returns the face area vector in cell i of face j" 
if (i < 0) or (i >= self .ncells) : 

print "Error: argument i out of range in geometry . avect" 
exit () 

if (j < 0) or (j >= self .nf aces) : 

print "Error: argument j out of range in geometry .avect" 

exitO 
ii = i*self .nf aces + j 
return self . area_vector [ii ] 

# 

def adjvrtx(self , i, j , k) : 

■Returns the adjacent cell and vertex associated with cell i, vertex j, 

face k" 

if (i < 0) or (i >= self .ncells) : 

print "Error: argument i out of range in geometry .ad jvrtx" 
exitO 

if (j < 0) or (j >= self.nvrtx): 

print "Error: argument j out of range in geometry . ad jvrtx" 
exitO 

if (k < 0) or (k >= self .nf aces) : 

print "Error: argxmient k out of reuige in geometry .ad jvrtx" 
exit ( ) 

ii = i*self .nfaces*self .nvrtx + j*self.nvrtx + k 
return self .adjacent [ii] 

# 

def adjcell ( self , i , k) : 

"Returns the adjacent cell and face associated with cell i, face k" 
if (i < 0) or (i >= self . ncells) : 

print "Error: argument i out of range in geometry . ad j cell" 
exit < ) 

if (k < 0) or (k >= self .nf aces) : 

print "Error: argxament k out of range in geometry . ad j cell" 
exitO 

j = self . face_node[k] [0] 

ii = i*self . nf aces* self .nvrtx + j* self. nvrtx + k 
return ( self . adjacent [ii] [0] , self .adjacent [ii] [2] ) 

# 

def invariant ( self ) : 

"Checks the mesh data for internal consistency" 
if (len(self . cell_region) != self . ncells ) : 

print "Error: region data size != ncells" 
exit () 

if (len(self . cell_node) != self . ncells) : 

print "Error: celljiode data size != ncells" 
exit ( ) 

if (len(self .cell_vol\jme) != self .ncells) : 

print "Error: cell volume data size != ncells" 
exitO 

if (len(self .node_coord) != self .nnodes) : 

print "Error: node_coord data size != nnodes" 
exit ( ) 

if {len(self . face_area) != self . ncells*self .nfaces) : 
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print "Error: face area data size •= ncells*nf aces" 

exitO 

if ( len( self . area_vec tor ) != self . ncells*self .nf aces) : 

print "Error: face area vector data size != ncells*nf aces " 
exit () 

if (len(self .node_coord) != self .nnodes) : 

print "Error: node_coord data size != nnodes" 
exit ( ) 

if (len( self .adjacent) != self .ncells*self ,nvrtx*self . nf aces) : 
print "Error: adjacent pointer data size incorrect!" 
exit ( ) 

if ( min(self .cell^voliime) <= 0. ): 

print "Error: Zero or negative cell volumes found." 
exitO 

if < inin(self . face_area) <=0. ): 

print "Error: Zero or negative cell face areas found." 
exit ( ) 

for icell in xrange (self . ncells) : 

for ivrtx in xrange (self . nvrtx) : 

for iface in xrange (self . nf aces ) : 

XX = self . ad jvrtx( icell, ivrtx, iface) 
jcell = xx[0] 
jvrtx = XX [1] 
if (jcell < 0) : 

continue 
else: 

found = 0 

for i in xrange ( sel f .nf aces) : 

yy = self .adjvrtx( jcell, jvrtx, i) 
if (yy[03 == icell) and (yy[13 == ivrtx): 
jface = i 
found = 1 
break 
if ( found != 1 ) : 

print "Error: no matching pointer in adjacent array." 
exit ( ) 

if ( self . cell_node [jcell] [jvrtx] != 
self .cell_node[ icell] [ivrtx] ) : 
print "Error: cell node mis-matcli in adjacent array. " 
exit ( ) 

if ( self .adjcelKicell, iface) != { jcell, jface ) ): 
print "Error: adjacent cell pointers do not match" 
exit ( ) 

if ( self .adj cell (jcell, jface) != ( icell, iface ) ): 
print "Error: adjacent cell pointers do not match" 
exit { ) 
for i in xrange (3) : 

zz = self .avect (icell, iface) [i] + \ 

self .avect (jcell, jface) [i] 
if ( abs(zz) > l.e-15) : 

print "Error found in face area vectors" 
exit ( ) 

return 1 

# 

# end of geom.py 

# 

# 

# 

# File: matp.py 



# 

# Provides multi-group material properties 

# information. Anisotropic up, down, euid self scatter 

# are provided for. Anisotropic scattering 

# data is provided in terms of Legendre (Pn) moments of 

# the scattering cross sections. 
# 

# All groups cind materials are assumed to have the same 
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# Pn scattering order for this demo code. However, considerable 

# memory and computational expense may be saved in a production 

# code by allowing the Pn scattering order to vary from one 

# group to another within a material. Likewise, each material 

# could be allowed to have its own individual Pn order. 
# 



# Author: John McGhee, Radion Technologies 
# 

# Date : 19 Feb 2004 

# 

# 

# $Id: matp.py,v 1.6 2004/03/05 21:51:31 mcghee Exp $ 
# 

# 

from sys import exit 

class Materia 1 Props : 

"Provides material property data." 

def init (self, matprop_f ile) : 

"Reads material properties data from disk." 

from string import split, atof, atoi, strip 

# Setup for read of fixed matprop data from disk, 
self . name=matprop_f ile 

f = open (matprop_f ile) 
record = f.readlineO 

# Read dimensions of the material properties from disk, 
self .matprop = [] 

while (record and strip (record) != "dimensions" ): 

record = f . readline ( ) 
if (record == " " ) : 

print "Error: imable to find dimensions keyword on matprop file %s\n" % 

matprop_f ile 

exit ( ) 
record = f.readlineO 
j = split (record) 
if (len( j) != 4) : 

print "Error: Incorrect number of dimensions in matprop file" 
self .ngroups=atoi ( j [0] ) 
self.nscxs = atoi(j[l]) 
self.nmat = atoi{j[2]) 
self.icp = atoi(j[3]) 
record = f . readline ( ) 

while (record and strip (record) != " end_dimensions " ): 

record = f . readline ( ) 
if (record == " " ) : 

print "Error: unable to find end^dimensions keyword on matprop file 
%s\n" % matprop_file 
exit ( ) 

if (self .ngroups <= 0) or (self.nscxs <= 0) or (self.nmat <= 0): 
print "Error: dimensions out of range on matprop file %s\n" % 

matprop_f i le 

exitO 

# Read in the matprop values from disk. 

while (record and strip(record) != "mater ial_data" ): 

record = f . readline ( ) 
if (record == ""): 

print "Error: unable to find matprop keyword on matprop file %s\n" % 

matprop_f i le 

exit () 



self .matprop = [] 
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record = f . readline ( ) 

while (record and strip (record) != "end_material_data'' ) : 
j = split (record) 
if (len(j) == 0) : 

print "Error: blank line found in matprop data." 
exitO 

self .matprop. append (map (atof, j ) ) 
record = f . readline ( ) 
if (record == " " ) : 

print "Error: unable to find end_matprop keyword on matprop file %s\n'' 
% matprop_file 

exitO 

self .nmatprop = len (self .matprop) 

ii = (self . nscxs+1) *self .nmat*self .ngroups 
if (self .nmatprop != ii) : 

print "Error: incorrect number of matprop values on file" 

# 

def sigsct(self, 1, isrc, isink, imat) : 
"Gets scattering cross section." 
if (imat < 0) or (imat >= self .nmat) : 

print "Error: imat argiiment out of range in scat" 
exit () 

if (isink < 0) or (isink >= self .ngroups) : 

print "Error: isink arg\ament out of range in scat" 
exit ( ) 

if (isrc < 0) or (isrc >= sel f . ngroups ) : 

print "Error: isrc argument out of range in scat" 

exit () 
if (1 < 0) : 

print "Error: 1 argument out of range in scat" 

exit () 
if (1 >= self.nscxs): 

return 0 . 
else : 

ii = imat* self .ngroups* (self .nscxs+1) + isrc* (self .nscxs+1) +1+1 
return self .matprop [ii] [isink] 

# 

def sigt(self , igroup, imat) : 

■Gets total cross section. " 

if (imat < 0) or (imat >= self .nmat) : 

print "Error: imat argument out of range in sigt" 

exit ( ) 

if (igroup < 0) or (igroup >= self .ngroups ) : 

print "Error: igroup argxament out of range in sigt" 
exit () 

ii = imat* self .ngroups* (self .nscxs+1) + igroup* ( self . nscxs+1) 
return self .matprop [ii] [0] 



def sigreac (self , igroup, imat): 

"Gets reaction cross section, " 

if (imat < 0) or (imat >= self .nmat) : 

print "Error: imat argument out of range in reac" 
exit ( ) 

if (igroup < 0) or (igroup >= self .ngroups) : 

print "Error: igroup argument out of range in reac" 
exitO 

ii = imat*self .ngroups* (self .nscxs+1) + igroup* (self .nscxs+1) 
return self .matprop [ii] [1] 



# 
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def sigstplself , igroup/ imat) : 

"Gets the stopping power over delta E." 
if (imat < 0) or (imat >= self .nmat) : 

print "Error: imat argument out of range in sigstp" 
exit ( ) 

if (igroup < 0) or (igroup >= self . ngroups) : 

print "Error: igroups argxjment out of range in sigstp" 
exitO 

ii = iinat*self .ngroups* (self .nscxs+1) + igroup* (self .nscxs+1) 
return self .inatprop[ii] [2] 

# 

# end of matp.py 

# 

# Multi-group Material Properties data 
# 

# Format : 

# units are 1/cm for cross sections, 

# MeV for energy, scattering cross sections 

# are given in Legendre moments. 
# 

# dimensions ngroups, nscxs, nmaterials, icp 
# 

# Material mat=0 

# Group g=0 

# sigma-t sig-reac sigma-s topping-power /delta-E 

# sigma-0,g->l, sigma-0 , g->2 , sigma-0 , g->3 . . . . sigma-0, g->ngroups 

# sigTna-l,g->l, 

# . 

# . 

# sigma-nscxs,g->l 
# 

# Group g=l 

# etc. 
# 

# Material mat=l 

# etc. 

dimensions 

2 2 2 0 
end_dimensions 

mater ial_data 

2.00 0.50 0.10 

0.50 0.20 

0.05 0.02 

4.00 1.00 0.20 

0.00 0.40 

0.00 0.04 

1.00 0.25 0.05 

0.25 0.10 

0-025 0.01 

2.00 0.50 0.10 

0.00 0.20 

0.00 0.02 
end_ma t er i a l_da t a 

# mesh geometry description 

# a (1x1x1) cube centered at (0.5, 0.5, 0.5) with 

# 2 regions, 862 tet cells, 207 nodes, and 200 triangular 

# boundary faces. Node numbering is zero based. Mesh regions 

# are assumed to be id'd with consecutive integers beginning 

# with 0. 

nodes 

0 0 . OOOOOOOOOOE+00 4 . OOOOOOOOOOE-01 0 . OOOOOOOOOOE+00 

1 0. OOOOOOOOOOE+00 1. OOOOOOOOOOE+00 1. OOOOOOOOOOE+00 
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2 l.OOOOOOOOOOE+00 

3 l.OOOOOOOOOOE+00 

4 O.OOOOOOOOOOE+00 

5 O.OOOOOOOOOOE+00 

6 l.OOOOOOOOOOE+00 

7 l.OOOOOOOOOOE+00 

8 O.OOOOOOOOOOE+00 

9 4.0000000000E-01 

10 4.0000000000E-01 

11 O.OOOOOOOOOOE+00 

12 4.0000000000E-01 

13 4.0000000000E-01 

14 O.OOOOOOOOOOE+00 

15 O.OOOOOOOOOOE+00 

16 2 .OOOOOOOOOOE-01 

17 2 .OOOOOOOOOOE-Ol 

18 O.OOOOOOOOOOE+00 

19 2. OOOOOOOOOOE-Ol 

20 2. OOOOOOOOOOE-Ol 

21 O.OOOOOOOOOOE+OO 

22 O.OOOOOOOOOOE+OO 

23 l.OOOOOOOOOOE+00 

24 l.OOOOOOOOOOE+00 

25 l.OOOOOOOOOOE+00 

26 l.OOOOOOOOOOE+00 

27 l.OOOOOOOOOOE+OO 

28 l.OOOOOOOOOOE+OO 

29 O.OOOOOOOOOOE+OO 

30 O.OOOOOOOOOOE+OO 

31 O.OOOOOOOOOOE+OO 

32 6.8750000000E-01 

33 3 .8671875000E-01 

34 1.9335937500E-01 

35 l.OOOOOOOOOOE+OO 

36 l.OOOOOOOOOOE+OO 

37 l.OOOOOOOOOOE+OO 

38 3 .1250000000E-01 

39 5.5389985616E-01 

40 8. 0664062500E-01 

41 O.OOOOOOOOOOE+OO 

42 O.OOOOOOOOOOE+00 

43 O.OOOOOOOOOOE+OO 

44 6.8750000000E-01 

45 3.8671875000E-01 

46 1.9335937500E-01 

47 l.OOOOOOOOOOE+OO 

48 1. OOOOOOOOOOE+OO 

49 l.OOOOOOOOOOE+OO 

50 7. OOOOOOOOOOE-Ol 

51 O.OOOOOOOOOOE+OO 

52 4. OOOOOOOOOOE-Ol 

53 4. OOOOOOOOOOE-Ol 

54 O.OOOOOOOOOOE+00 

55 2. OOOOOOOOOOE-Ol 

56 4. OOOOOOOOOOE-Ol 

57 2. OOOOOOOOOOE-Ol 

58 O.OOOOOOOOOOE+OO 

59 2. OOOOOOOOOOE-Ol 

60 4. OOOOOOOOOOE-Ol 

61 2. OOOOOOOOOOE-Ol 

62 2. OOOOOOOOOOE-Ol 

63 4. OOOOOOOOOOE-Ol 

64 O.OOOOOOOOOOE+00 

65 O.OOOOOOOOOOE+OO 

66 O.OOOOOOOOOOE+00 

67 0. OOOOOOOOOOE+OO 

68 O.OOOOOOOOOOE+00 



l.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
4. OOOOOOOOOOE-Ol 
4. OOOOOOOOOOE-Ol 
O.OOOOOOOOOOE+OO 

O.OOOOOOOOOOE+OO 
4. OOOOOOOOOOE-Ol 
O.OOOOOOOOOOE+OO 
2. OOOOOOOOOOE-Ol 
2. OOOOOOOOOOE-Ol 
2. OOOOOOOOOOE-Ol 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
2. OOOOOOOOOOE-Ol 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 



l.OOOOOOOOOOE+OO 
3.1250000000E-01 
6.1328125000E-01 
8.0664062500E-01 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
6.8750000000E-01 
1582841730E-01 
9335937500E-01 
OOOOOOOOOOE+OO 
OOOOOOOOOOE+OO 
OOOOOOOOOOE+OO 
1250000000E-01 
1328125000E-01 
8.0664062500E-01 
O.OOOOOOOOOOE+OO 



.OOOOOOOOOOE-Ol 
. OOOOOOOOOOE+OO 
. OOOOOOOOOOE-Ol 
.OOOOOOOOOOE-Ol 
. OOOOOOOOOOE-Ol 
2. OOOOOOOOOOE-Ol 
O.OOOOOOOOOOE+00 
.OOOOOOOOOOE-Ol 
.OOOOOOOOOOE-Ol 
.OOOOOOOOOOE-Ol 
.OOOOOOOOOOE-Ol 
.OOOOOOOOOOE-Ol 
.OOOOOOOOOOE-Ol 
.9816277408E-01 
. 6429795369E-01 
6.9020637049E-01 
6.8143850763E-01 
8.5642346013E-01 



l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+OO 
0. OOOOOOOOOOE+OO 
4. OOOOOOOOOOE-Ol 
4. OOOOOOOOOOE-Ol 
4. OOOOOOOOOOE-Ol 

4. OOOOOOOOOOE-Ol 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
2. OOOOOOOOOOE-Ol 
2. OOOOOOOOOOE-Ol 
2. OOOOOOOOOOE-Ol 
2. OOOOOOOOOOE-Ol 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
7. OOOOOOOOOOE-Ol 
6.8750000000E-01 
3.1175142799E-01 
1.9335937500E-01 
6.8750000000E-01 
3.8671875000E-01 
1.9335937500E-01 
6.8750000000E-01 
3.8671875000E-01 
1.9335937500E-01 
l.OOOOOOOOOOE+00 
1. OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+00 
1. OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+00 
1. OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
2. OOOOOOOOOOE-Ol 
2. OOOOOOOOOOE-Ol 
2. OOOOOOOOOOE-Ol 
4. OOOOOOOOOOE-Ol 
4. OOOOOOOOOOE-Ol 
4. OOOOOOOOOOE-Ol 
4. OOOOOOOOOOE-Ol 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
4. OOOOOOOOOOE-Ol 
2. OOOOOOOOOOE-Ol 
2. OOOOOOOOOOE-Ol 
8.5265650289E-01 
6.0022303814E-01 
8.1647869761E-01 
5.0305088449E-01 
2.0469430523E-01 
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O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+00 
5.4286513258E-01 
1.6114548438E-01 
3.0058947384E-01 
3.8506624883E-01 
7 .9513906813E-01 
6.3890496821E-01 
8.7377401275E-01 
8.5520505408E-01 
5.1427246296E-01 
2.0436734847E-01 
2.1056348860E-01 
5.4576935545E-01 
8.1762242372E-'01 
7.5109655258E-01 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
1.7294999713E-01 
8.5207743342E-01 
2.0649949350E-01 
8.0865970785E-01 
5.0279034487E-01 
2.5161078888E-01 
6.8115234375E-01 
5.3393554688E-01 
9 .3139665146E-02 
1.8632077853E-01 
2.0524985322E-01 
8 . 0539233807E-01 
8.5732329126E-01 
2.3515465154E-01 
4.8698290591E-01 
7.6836884278E-01 
5.3945232577E-01 
5.3129423908E-01 
5.2142553338E-01 
6.9057450500E-01 
6.3589775149E-01 
5.2610324720E-01 
5.3892400556E-01 
6.6999320000E-01 
6.4736087761E-01 
1. 6280610078E-01 
1.5568563329E-01 
3 .3071891700E-01 
2.8757807380E-01 
1. 6153279295E-01 
1.7051505099E-01 
3 .4376014959E-01 
2.9479372820E-01 
5,1698251895E-01 
6.7969530916E-01 
5.1618508017E-01 
6.8977402522E-01 
4.5891086203E-'01 
4.6699594988E-01 
6.6482681022E-01 
6.5235353985E-01 
1.5805075826E-01 
1.6877853480E-01 
3.2986998766E-01 



5.4274727097E-01 
1.5658658666E-01 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+OO 
7.9331699083E-01 
6.9535867584E-01 
8.5696128358E-01 
5.4411661671E-01 
2.0873396013E-01 
1.8451589823E-01 
4.9052150087E-01 
8.2402927123E-01 
1.2487433107E-'01 
7.9285776851E-01 
2.0144018136E-01 
4.8428942222E-01 
7.1967437916E-01 
2.7808336078E-01 
4.6606445313E-01 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+OO 
1.2553485179E-01 
2.6437242858E-01 
8.5230823230E-01 
1.4690864190E-01 
7.9149776183E-01 
5.2583210433E-01 
2 . 5022079020E-01 
4.4525994821E-01 
6.6623636466E-01 
3.6301742777E-01 
1.9223792132E-01 
1.7024114100E-01 
3.3465936049E-01 
3.5006957730E-01 
1.7231611860E-01 
1,5999830000E-01 
3.0573460412E-01 
7,3718449692E-01 
5.5970950946E-01 
5.5884902893E-01 
7 .0027992714E-01 
6.9803521734E-01 
5.5711183509E-01 
5.4232782442E-01 
6.6063039019E-01 
5.5944707707E-01 
5.2274344087E-01 
5.2586645378E-01 
4.7535559560E-01 
7.4247179192E-01 
6.8829280667E-01 
7.0275413419E-01 
6.4410824466E-01 
3 .5640783408E-01 
1.8716897893E-01 
1.7349008539E-01 



2.1012134137E-01 
5.4092971210E-01 
1.5566251981E-01 
5.3892617058E-01 
8.1720302964E-01 
5.8953247729E-01 
8.2792057671E-01 
3.8698326115E>01 
2.0445757272E-01 
O.OOOOOOOOOOE+00 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
O.OOOOOOOOOOE+OO 
1.7413106030E~01 
2.0700356803E-01 
8.5510723799E-01 
7.9444561904E-01 
2.5701066922E-01 
5.1294207125E-01 
4.7588364265E-01 
7.5459188612E-01 
1.7572284990E-01 
2.0504346565E-01 
8.5711728089E~01 
8.1686297897E-01 
2.4524971658E-01 
5.2225394214E-01 
5.1855468750E-01 
7.5927734375E-01 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+OO 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+00 
l.OOOOOOOOOOE+OO 
1.5615253846E-01 
1. 6863581000E-01 
1.9396915061E-01 
1.5092855295E-01 
3.20264e4946E-01 
3.2929723165E-01 
3.2999660000E-01 
3.2217274456E-01 
1.6134962498E-01 
1.5509874015E-01 
1.6928060623E-01 
1.5926615033E-01 
2.8470005310E-01 
3 .2806298508E-01 
3 .4041747705E-01 
2 .9101613 137E-01 
1.8200762808E-01 
1.8362101945E-01 
3.3424519004E-01 
3.0938541431E-01 
2.0331723849E-01 
2.9446984286E-01 
1.7337280910E-01 
3 .0403009171E-01 
5.2667402848E-01 
5.2041032434E-01 
5.3571406758E-01 
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137 


3 


.2196656960E- 


-01 


3 


.4774897588E- 


-01 


5 


.2691975944E- 


-01 


138 


1 


. 6241740361E- 


-01 


3 


. 1114277277E- 


-01 


6 


.4321023544E- 


-01 


139 


1 


.7351139147E- 


-01 


1 


.9251167902E- 


-01 


6 


,8578253762E- 


-01 


140 


3 


,5033875261E- 


-01 


1 


. 6496039409E- 


-01 


6 


.8906200411E- 


-01 


141 


3 


.2397395872E- 


-01 


3 


.0700405926E- 


-01 


6 


.4879962237E- 


-01 


142 


5 


.0425714015E- 


-01 


3 


.3498460286E- 


-01 


5 


.0446401895E- 


-01 


143 


5 


.2240052366E- 


-01 


1 


.7979297462E- 


-01 


5 


.2044836141E- 


-01 


144 


7 


.0030733512E- 


-01 


1 


.9580382933E- 


-01 


5 


.3858694926E- 


-01 


145 


6 


.8396487320E- 


-01 


3 


.1011976595E- 


-01 


5.1460594801E- 


-01 


146 


5 


.2994296233E- 


-01 


1 


.6883896656E- 


-01 


7 


3688375136E- 


-01 


147 


5 


.1980765072E- 


-01 


3 


.0646663921E- 


-01 


6 


8667479623E- 


-01 


148 


6 


.8799537735E- 


-01 
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This a multi-group 3D arbitrary tetrahedral mesh discrete ordinates 
demonstration code written in Python. The purpose of this code is to 
demonstrate the general features of the solution algorithm. It is not 
5 optimized for computational efficiency and is not intended to be used 
for large problems or in a production environment. 

Written by: 

10 John McGhee, Radion Technologies, Los Alamos NM, 19 Feb 2004. 

This program is written in Python. Python is a widely availcUDle, 
powerful, easy to use, and easy to learn • scripting language. 
(And it ' s free ! ) References : 
15 "Learning Python", Mark Lutz and David Ascher, O'Reilly, 1999. 

— "Python Essential Reference", David Beazley, New Riders, 2000. 
Python Home Page — http: //http: / /www. python. org/ 
The Vaults of Parnassus — http://www.vex.net/parnassus/ 
-- Numerical Python -- http://www.pfdubois.com/numpy/ 
20 -- Scipy -- http://www.scipy.org/ 

Starship Python - http://starship.python.net 

A discussion of the design philosophy of the code and various other 
background information is provided in the file background_inf o . asc . 
25 Extensive comments are provided in the body of the source code that 
explain details of the algorithm. 

A small sample problem is included for demonstration purposes. All 
output is written to the screen with the exception of a GMV 
30 visualization link file. Sample screen output is included in the file 
"frost. out_std" . A coxi^ressed GMV output file is included as 
"gmv. out_std.gz" . GMV is a publicly available visualization utility 
that can be obtained from Los Alamos National Laboratory, Los Alamos NM. 
See http://www-xdiv.lcUil.gov/XCM/gmv/ 

How to run the code: 

Simply type . /frost. py at the command prompt to execute the included 
sample problem. Runs in about two minutes, sending output to the screen as 
the calculation progresses. Other problems can be run by replacing the 
40 data files: 

matprop . inp , 

aquad . inp , and , 
mesh . inp 

45 with the data files of choice and redefining the fixed source as desired 
in the file "fsrc.py". Formats for the data files are briefly defined 
in the included demo files. Adding the command line option "fsds" will 
cause the first scattered distributed source solution demonstration mode 
to be invoked. 

50 

Limitations and things Frost will not do: 

tets only, no other spatial element types 
-- standard scattering only, no Galerkin scattering 
55 -- PI scattering only, no higher order scattering terms 
--no meshes which are not directed acyclic graphs (DAG) 

only vacuum boundary conditions 
-- no dsa of scattering source 

fsds option assumes a single isotropic point source. Only first 
60 scattered source problems are supported, no second or subsequent 

collided sources. 
-- no upscatter in energy, down-scatter only. 

only a one-to-one correspondence between mesh regions and materials. 

is supported, no material mixing. 

65 

Things that are included in Frost: 

-- Geometry Class - geom.py - Reads in mesh off disk, computes tet 
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volxjmes, face areas, face normals, connectivity, pointers to adjacent 
cells. Provides accessors methods to get to it all. 

— Angular Quadrature Class - aquad.py - Reads in quadrature off disk. 
Accessor methods to get at data. Computes discrete to moment and moment 
to discrete matrices. Only does standard scattering, no Galerkin. 

— Fixed Source Class - fsrc.py - Defines volumetric fixed sources. 
Provides accessor methods to get at the data. 

— BteEquation Class - bte.py - Forms left and right hand 

side of LD Boltzmann transport equation for neutral and charged 
particles. Solves the system with LU decomp. 

Frost - frost. py - Main driver. Calls all the methods described above, 
loops through groups, angles, and mesh cells. Computes baleuice table for 
6dit and debug. Writes visualization dump. 

-- Sweep Ordering Class - sord.py - Orders the mesh cells from upwind 
to downwind for any specified quadrature direction. Does not do graph 
cycles . 

— Spherical Harmonics methods - sph.py - Evaluates the orthonormal 
spherical harmonics surface fxinctions for any given direction. Only 
good through degree 1 . 

— Material Properties Class - matp.py - Reads in multi-group material 
properties from disk. Provides accessor functions to get at the data. 

— Edits Class - edit.py - Writes GMV link file for visualization. Useful 
for debug and other post-processing. 

— Scattering Methods - scat.py - Does out of group scatter and self scatter 
separately, iterates to convergence with good particle balance. 

— Charged Particles - This was accomplished as an addition 
to the BteEquation class. 

First Scattered Distributed Source Class - fsds.py provides an 
analytic method to calculate the uncollided solution component and 
computes a first scattered distributed source. Also provides a method 
to add the uncollided and collided components together to form the 
total solution. 

— Ray Tracing Class - trace. py provides algorithms for tracing rays 
on tetrahedral meshes, related algorithms for interpolation of 
fields, euid methods for finding cells which contain arbitrarily 
located spatial points. 



— jmm, Los Alamos, 10 Mar 2004. 

# 

# 

# File: scat.py 
# 

# Creates scattering source moments. 
# 

# Author: John McGhee, Radion Technologies 
# 

# Date : 19 Feb 2004 

# 

# 

# $Id: scat.py,v 1.6 2004/03/05 23:01:13 mcghee Exp $ 
# 

# 

from sys import exit 
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def in_scatter (icell, ig, meshdata, matprop, momdata, qdata) : 

"Produces the in_scatter source for cell icell and group ig) " 

result = [] 
isink = ig 

imat = meshdata. region (icell) 
for imom in range (qdata .nmom) : 
1 = qdata . index_lm [ imom] [ 0 ] 
X = [ 0., 0., 0., 0.] 
for isrc in range (matprop. ngroups) : 
if (isrc == isink) : 

continue 
for j in reuige ( 4 ) : 

x[j] = x[j] + momdata [imom] [isrc] [icell] [j ] * \ 
matprop. sigsct (1, isrc, isin]c, imat) 
resul t . append ( x ) 
return result 

# 

def self_scatter ( icell , ig, meshdata, matprop, momdata, qdata): 

"Produces the self_scatter source for cell icell eund group ig) " 

result = [] 
isrc = ig 
isinJc = ig 

imat = meshdata. region (icell) 
for imom in range (qdata. nmom) : 

1 = qdata , index_lm [ imom ] [ 0 ] 

X = [] 

for j in range (4) : 

X . append (momdata [ imom] [ ig ] [ i cell ] [ j ] * \ 

matprop. sigsctd/ isrc, isin]c, imat) ) 
result . append (x) 
return result 

# 

# end of scat.py 

# 

# 

# 

# File: sord.py 



# 

# Determines the ordering of the cell solution process 

# for a particular direction eind mesh. 
# 

# This demo assumes a tet mesh, but other arbitrary finite 

# element meshes can be ordered in a similar manner. For elements 

# with curvilinear faces, a face average normal is often used 

# to determine in/out flow status of the face. 
# 

# The Sn sweep dependencies can be defined in terms of a directed graph. 

# We define the vertexes of the graph to be the mesh cells, and the 

# edges of the graph to be the outgoing faces on a cell directed to 

# the incoming faces on the adjacent cells. We get a different graph for 

# every Sn eingle on every mesh. 
# 

# The process of finding a sweep ordering is what is Icnown as a 

# topological sort in graph theory. We hope that the sweep- ordering for 

# each angle can be represented as a directed-acyclic-graph (dag) . 

# If strongly connected components (cycles) are present, 

# then the mesh graph is not a dag and no topological sort is 

# possible. 
# 

# If cycles are discovered, a graph theory method Icnown as a 

# "depth- first-search" can be used to convert the graph into a 



Docket No. 200313554-1 



112 



# dag. A depth first search partitions the 

# graph edges into four types: tree-edges, f orward- edges , cross-edges, 

# and back-edges. Viewed from any vertex, the depth first search for 

# a dag has no back-edges. Thus the graph can be converted to a 

# dag by removing the back edges. 
# 

# For purposes of simplification, the algorithm herein assumes 

# that the mesh graph is a dag. 
# 

# 

# Author: John McGhee, Radion Technologies 
# 

# Date : 19 Feb 2004 

# 

# 

# $Id: sord.py,v 1.6 2004/03/05 21:51:31 mcghee Exp $ 
# 

# 

from sys import exit 



class SolutionOrder : 

"Provides information on the mesh cell ordering required for solution," 

def init (self, geom, cjuad) : 

"Sets up the sweep ordering" 
self.ncells = geom.ncells 
self.neuig = guad.nang 
self. sweep_or der = [ ] 

# Loop over angles 

for kk in xrange (quad. nang) : 

nnd = self.ncells 

omega = quad . qdpt ( kk ) 

self . sweep_order . append ( [ ] ) 

f aces_needed = [ ] 

out_face = [] 

# Find the outgoing cell faces for this angle, 
for icell in xrange ( self . ncells) : 

i = 0 

out_f ace . append ( [ ] ) 
for iface in range (geom. nf aces) : 
XX = 0. 

yy = geom. avect (icell, iface) 
for j in remge (geom.ndim) : 

XX = XX + omega[ j ] *yy [j] 
if ( ( XX < -geom.dplimit ) and 

{ geom. adjcell (icell, if ace) [0] != geom. bdry_f lag [0] ) ): 
i = i+1 
elif ( XX > geom.dplimit ) : 

out_f ace [icell] . append (if ace) 
f aces__needed . append ( i ) 
if { i == 0 ) : 

self . sweep_order [kk] . append ( icell ) 
nnd = nnd - 1 

# Go through the mesh in a wavefront manner, ordering the 

# cells by wave front number, 
i = 0 

this__wave = self . sweep_order [kk] [ : ] 
if (len{this_wave) == 0) : 

print "Error: uneible to find any starting cells." 

print "kk= ", kk, ", nnd= ", nnd, ", swp= i 

exit ( ) 
while (nnd > 0) : 
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i = i+1 

if ( i >= self.ncells) : 

print "Error: niimber of sweeps exceeds ncells!" 
exit ( ) 
next_wave = [ ] 
for icell in this_wave: 

for iface in out_f ace [icell] : 

jcell = geom. adj cell ( icell , if ace) [0] 
if (jcell == geom. bdry_f lag [ 0] ) : 
continue 

faces_needed[ jcell] = faces_needed[ jcell] - 1 

if <faces_needed[ jcell] == 0) : 
next_wave . append ( j cell ) 
self . sweep_order [kk] . append ( j cell ) 
nnd = nnd - 1 
if (len (next_wave) == 0) : 

# At this point a cycle in the mesh graph has been 

# found. A production algorithm would call on a 

# method to break this cycle by flagging back edges 

# in the mesh for iteration. 

print "Error: graph cycle found, unable to continue sweep 

ordering process . " 

print "kk= " , kk, nnd= nnd, " , swp= " , i 
exit ( ) 
this_wave = next_wave 

if ( len(self . sweep_order [kk] ) != self.ncells ): 
print "Error: did not order all cells." 
exit ( ) 

# return for next angle. 

# 

def order (self , kk) : 

"Returns the cell-wise order of solution for angle kk" 
if (kk < 0) or (kk > self .nang) : 

print "Error: angle argument out of range in cell^order" 
exit ( ) 

return self . sweep_order [kk] 

# 

# end of sord.py 

# 

# 

# 

# File: sph.py 



# 

# Evaluates the orthonormal spherical harmonics surface 

# functions. This function is defined in standard 

# mathematical references. Normalization used herein 

# assumes integration over the unit sphere = 1. 

# This demo function generator only supports degree 

# <= 1, a production version would support any 

# requested degree and order. 
# 

# 1 is the degree of the function (1>=0) 

# m is the order of the function (-l<=m<=l) 

# xmu is the cosine of the polar angle. The azimuthal angle 

# phi is measured in the plus xi direction from the plus eta axis. 

# sph = [p(l,abs(m) ,xmu) ]*[cos(abs(m)*phi) ] for m.ge.O; 

# = [p(l,abs (m) ,xmu) ] * [sin(eJDS (m) *phi) ] form. It. 0 

# where p(l,m,xmu) is the Associated Legendre Polynomial 

# m = 0 returns the Legendre Polynomials 
# 
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# Author: John McGhee, Radion Technologies 
# 

# Date : 19 Feb 2004 

# 

# 

# $Id; sph.py,v 1.6 2004/03/05 21:51:31 mcghee Exp $ 
# 

# 

def sphf\jn(l, m, xmu, eta, xi) : 

"Evaluates the ortho-normal spherical harmonics surface fxinction" 

from sys import exit 

from math import sqrt, cos, sin, acos, pi 

if ( 1 < 0 ) : 

print "Error; 1 must be >= 0 in sph" 

exit ( ) 
if ( abs(in) > 1 ) : 

print "Error: abs(m) must be <= 1" 

exit ( ) 

if (1 == 0) : 

return 1 . 
elif (1 == 1) : 

if (m == 0) : 

return sgrt ( 3 . ) *xmu 

else: 

# Catch the limiting case of xmu = +1,-1 

XX = min ( ( 1 . -xmu) * ( 1 . +xmu) , (eta*eta + xi*xi) ) 
if (XX < l.e-16) : 
return 0 . 

# Get the azimuthal angle, phi 
XX = eta*eta 

XX = XX/ (xx + xi*xi) 

XX = max( 0., min(l., xx) ) 

XX = acos ( sqrt(xx) ) 

if (xi > 0.) : 

if (eta > 0. ) : 

pass 
else: 

XX = pi - XX 

else: 

if (eta > 0. ) : 

XX = 2 . *pi - XX 
else: 

XX = pi + XX 

# Return the function value 

yy = sqrt ( (1.- xmu)*(l.+ xmu) ) 
if (m == 1) : 

return sqrt (3 . ) *yy*cos (xx) 
else: 

return sqrt (3 , ) *yy*sin(xx) 

else : 

print "Error: this version of sph only supports 1=0,1" 
exit () 



# 
# 
# 



end of sph.py 



